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A  solution  of  groups  technique  was  developed  for  use 
with  fluctuation  solution  theory.   The  general  expressions 
for  calculation  of  pressure  and  chemical  potential  changes 
from  some  fixed  reference  states  have  been  shown.   A  new 
corresponding  states  theory  correlation  for  direct  correla- 
tion function  integrals  was  proposed  and  used  with  the 
group  contribution  technique  for  calculation  of  pressure 
changes  during  compression  for  several  n-alkanes  and 
methanol. 

This  work  gives  a  detailed  analysis  of  the  RISM  theory 
of  liquids.   Shown  are  new  results  for  perturbation  theory 
and  a  generalized  compressibility  theorem  for  RISM  fluids. 
The  use  of  the  RISM  theory  for  calculation  of  thermodynamic 
properties  of  real  fluids  also  is  given. 


The  use  of  hard  sphere  reference  fluids  for  development 
of  equations  of  state  has  been  explored.   A  generalized 
hard  sphere  equation  of  state  was  developed.   It  was  shown 
that  the  most  accurate  hard  sphere  equation  of  state  is 
not  the  best  reference  system  for  construction  of  a  liquid 
phase  equation  of  state  of  the  van  der  Waals  form. 


CHAPTER  1 
INTRODUCTION 


The  rational  design  of  chemical  process  equipment 
requires  knowledge  of  the  thermodynamic  properties  of  the 
substances  involved.   More  specifically,  volumetric 
properties  are  needed  to  size  a  piece  of  equipment, 
enthalpies  are  needed  to  determine  heat  duties,  and  fugaci- 
ties  are  used  to  decide  feasibility  of  reactions  and 
separations . 

The  general  problem  in  physical  property  correlation 
and  estimation  is  then  to  determine  a  procedure  for  calcula- 
tion of  volumes,  enthalpies,  and  fugacities,  for  multi- 
component,  multiphase  mixtures.   This  has  been  accomplished 
for  some  pure  components  using  an  equation  of  state  with 
a  large  number  of  parameters.   Here,  a  more  modest  problem 
is  addressed.   A  formalism  is  developed  that  allows  the 
volumetric  properties  and  chemical  potentials  for  liquid 
mixtures  to  be  expressed  in  terms  of  a  single  group  of 
parameterized  functions. 

The  approach  used  here  tries  to  recognize  the  essential 
differences  between  liquid  and  vapor  phase  properties. 
A  thermodynamically  consistent  formulation  for  the  volumetric 
properties  and  chemical  potentials  is  developed  using  an 
equation  of  state.   Here,  unlike  many  other  schemes,  the 
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equation  of  state  uses  a  liquid  reference  state.   This 
is  all  made  possible  by  utilization  of  the  Kirkwood-Buf f 
(1951)  or  fluctuation  solution  theory. 

In  an  effort  to  allow  for  extrapolation  of  present 
results  to  state  conditions  not  examined,  the  correlations 
are  put  into  a  corresponding  states  form.   This  can  be 
justified  from  analysis  of  the  microscopic  theory  that 
is  the  basis  for  this  macroscopic  approach. 

Finally,  to  allow  for  actual  physical  property  estima- 
tion, the  theory  is  cast  in  a  "group-contribution"  form. 
This  is  easily  accomplished  formally  because  the  fluctuation 
solution  theory  relates  molecular  physical  properties  to 
microscopic  correlations.   These  quantities  are  well  defined 
for  groups,  or  sites,  as  they  will  be  often  termed. 

The  next  chapter  in  this  work  goes  into  detail  in 
discussing  the  thermodynamic  aspects  of  the  problem  and 
the  possible  approaches  that  are  available  to  express  the 
required  physical  properties.   Then,  the  flucutation  solution 
theory  is  outlined. 

Chapter  3  is  devoted  to  the  discussion  of  group  con- 
tribution modeling.   Two  possible  means  of  linking  the 
group  contribution  ideas  to  the  fluctuation  theory  are 
presented  and  their  differences  expounded  upon.   This  chapter 
concludes  with  the  formulation  of  the  physical  property 
relationships  from  group  correlation  function  integrals. 

Approaches  to  modeling  correlation  functions  are  the 
emphasis  of  Chapter  4.   General  forms,  as  suggested  by 


3 
microscopic  perturbation  theories,  are  outlined  and 
discussed.   The  thermodynamic  consistency  requirements 
for  model  formulation  are  given  along  with  the  final  expres-. 
sions  for  the  thermodynamic  property  changes  from  the  chosen 

form. 

Because  the  correlation  function  integrals  contain 
unknown  functions,  expressed  in  corresponding  states  form, 
some  experimental  data  are  required  for  the  model  parameteri- 
zation.  In  Chapter  5  possible  data  for  this  parameterization 
are  discussed  critically,  and  two  sets  of  general  correla- 
tions are  developed.   One  of  these  correlations  is  then 
decided  upon  based  on  the  problem  of  interest  and  accuracy. 

Chapter  6  then  presents  calculations  for  several 
n-alkane  molecules'  compressions.   Both  correlations  and 
predictions  are  shown.   Chapter  7  is  a  general  discussion 

of  all  of  the  previously  reported  work,  and  Chapter  8  offers 

several  conclusions  and  recommendations  for  future  work. 
There  are  also  several  appendices  which  give  greater 

detail  on  developments  presented  in  the  body  and  listings 

of  several  useful  computer  programs. 


CHAPTER  2 
THERMODYNAMIC  THEORY 

Thermodynamics  is  a  science  that  has  its  basis  on 
certain  empirical  observations  that  have  become  known  as 
laws.   The  simplest  of  these  observations,  and  the  ones 
most  often  used  in  practice,  are  that  mass  and  energy  are 
conserved  quantities.   These  statements  allow  for  the  con- 
struction of  useful  mass  and  energy  balance  equations. 
These  balances,  when  augmented  with  constitutive  relations 
such  as  the  conditions  of  phase  or  reaction  equilibria, 
can  be  used  for  design  of  chemical  process  units.   However, 
the  power  of  these  balance  expressions  cannot  be  fully 
utilized  unless  the  physical  properties  that  appear  in 
the  expressions  are  available. 

In  the  following  discussion  one  other  empirical  obser- 
vation is  required.   This  is  not  known  as  a  law  but  is 
often  listed  as  a  postulate  of  thermodynamics  (Modell  and 
Reid,  1974).   This  result  may  be  stated  in  many  ways,  but 
the  simplest  formulation  is  as  follows:   For  a  simple, 
homogeneous  system  of  N  components,  with  only  P-V  work 
and  thermal  interactions  with  its  surroundings,  that  N+l 
independent,  intensive  variables  are  required  to  specify 
the  state  of  the  system.   The  simplest  use  of  this  observa- 
tion is  the  existence  of  equations  of  state,  such  as 
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P  =  f(T,V,X)  '        (2-1) 

Thermodynamic  Properties  of  Interest 

As  stated  previously,  this  work  is  concerned  with 
the  volumetric  properties  and  chemical  potentials  of  the 
components  in  liquid  mixtures.   The  chemical  potentials 
are  of  interest  for  phase  or  chemical  equilibrium  calcula- 
tions because  of  the  constraints  that  are  imposed.   For 
equilibrium  between  phases  a  and  6  the  constraints  are 

ya  =  uB  V.  in  both  a  and  3  (2-2) 

1     ix 

and  the  chemical  equilibrium  constraints  are 

Y  v. .  u.  =  0    V   reaction  j  (2-3) 

l 

These  equilibrium  conditions  are  used  to  determine  both 
feasible  operating  conditions  for  chemical  processes  and 
minimum  work  requirements  for  some  processes. 

To  size  a  piece  of  process  equipment  one  must  have 
knowledge  of  volumetric  properties  of  the  fluids  involved, 
The  molar  volume  of  a  mixture  is  found  from  knowledge  of 
the  component  partial  molar  volumes  by 


V  =  1  x.  V, 

l 


(2-4) 
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An  important  aspect  that  must  be  considered  in  con- 
struction of  thermodynamic  models  is  the  relationship  between 
chemical  potentials  and  partial  molar  volumes 

(  !H  )     =  v.  (2-5) 

3P    T,X     X 

When  both  the  chemical  potentials  and  partial  molar  volumes 
are  determined  from  the  same  equation  of  state,  then  equation 
(2-5)  will  always  be  satisfied.   However,  if  different 
correlations  are  used  for  calculating  these  two  different 
properties,  then  a  thermodynamic  inconsistency  exists  and 
can  cause  computational  difficulties. 

Calculation  of  Liquid  Volumes 

Reid,  Prausnitz,  and  Sherwood  (1977)  give  an  extensive 
review  of  techniques  for  correlating  liquid  molar  volumes. 
Most  effort  has  been  placed  on  calculation  of  saturated 
liquid  volumes  because  of  the  insensitivity  of  liquid  volumes 
to  pressure.   Most  of  the  correlations  are  in  corresponding 
states  form. 

Calculations  of  compressed  (subcooled)  liquid  volumes 
often  are  based  on  a  knowledge  of  the  saturated  liquid 
volume  at  the  system  temperature.   In  general,  the  volumes 
are  found  from  expressions  of  the  form 

V  =  V(T,P)  (2-6a) 


Vmix  =  V(T,P,X)  (2-6b) 

These  are  calculationally  convenient  forms.   If  a  complete 
equation  of  state  is  used,  the  volumes  are  found  from  solu- 
tion of  the  implicit  relation 

P  =  P(V,T,X) 

and  care  is  often  necessary  to  ensure  that  the  proper  volume 
is  calculated.   Equation  2-1  can  be  satisfied  for  several 
values  of  the  volume  at  a  fixed  temperature  and  composition. 
Partial  molar  volumes  are  found  from  the  definition 

v±  -  (  w:  )  (2-7> 

1   T,P,N  ... 

For  most  correlations  this  is  evaluated  by  using  one  of 
the  pure  component  correlations  with  a  set  of  mixing  rules 
applied  to  the  parameters. 

Chemical  Potentials  and  Fugacities 

Chemical  potentials  are  used  for  solution  of  phase 
or  reaction  equilibria  problems.   However,  the  chemical 
potentials  themselves  are  not  often  used  because  the 
equilibrium  expressions  can  be  written  in  terms  of  the 
fugacities,  defined  by 


U  =  y°  +  RT  in    (f\/f°)  (2-8) 


where  the  superscript  o  refers  to  a  given  reference  state 
of  a  specified  pressure,  composition,  and  phase  at  the 
system  temperature,  T.   The  fugacity  can  be  found  from 
two  general  approaches  based  on  different  reference  states. 
For  the  ideal  gas  reference  state  a  complete  equation  of 
state  is  required,  one  that  can  reasonably  predict  the 
system  volume.   Another,  more  common,  approach  is  to  use 
a  liquid  reference  state  and  an  expression  for  the  excess 
Gibbs  free  energy.   In  this  formulation  the  fugacity  is 
written  as 


.  =  X.  y.    f  . 
1     i'ii 


f  .  = 


2-9) 


The  Y.  term,  known  as  the  activity  coefficient,  is  found 
from  the  free  energy  model  and  is  used  to  correct  for  com- 
position nonideality  in  the  liquid  solution. 

Mathias  and  O'Connell  (1981)  have  proposed  a  slight 
variant  on  this  liquid  reference  state  scheme.   Their 
approach  is  based  on  using  the  temperature  and  the  component 
densities,  P.  =  N./V  as  the  independent  variables  to  describe 
the  state  of  the  system.   Then,  at  constant  temperature 
the  following  relation  holds  : 

3  Jin  f. 

Fugacity  ratios  can  be  calculated  based  on  any  reference 
state  if  the  partial  derivatives  in  equation  2-10  are  known. 


9 
Expressions  for  these  derivatives  in  terms  of  integrals 
of  microscopic  correlation  functions  are  presented  in  the 
next  section. 

Fluctuation  Solution  Theory 

Fluctuation  solution  theory  (Kirkwood  and  Buff,  1951; 
O'Connell,  1971,  1982)  is  a  bridge  that  connects  the  thermo- 
dynamic derivatives  to  statistical  mechanical  correlation 
function  integrals.   The  basic  relation  of  fluctuation 
theory  is 

3  <N.  > 
f  i  I  =  <N.N.>  -  <N.>  <N.>        (1    ,  ., 

where  the  brackets  denote  an  average  over  an  equilibrium 
grand  canonical  ensemble  and  B  =  lAgT.   These  averages 
are  related  to  correlation  function  integrals  by 


<N.N.>  -  5..  <N.>  =  C-K1)    J  dld2g.  .(1,2)     (2-12 

i  3  iD    i       o2  ^ 


where  Jdl  is  an  integration  over  all  phase  space  coordinates 
required  for  molecule  1  and  Q    is  the  normalization  constant 
for  the  orientation  dependence.   The  function  g^d/2) 
is  known  as  a  pair  correlation  function  and  is  directly 
related  to  a  two  molecule  conditional  probability  density. 
It  is  often  more  convenient  to  work  with  the  total  correla- 
tion function,  1^.(1,2),  defined  by 
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h.  .(1,2)  =  gij(l,2)  -  1  ■  (2-13) 


Combination  of  equations  2-12  and  2-13  leads  to 


3<N.> 

h—}  =    6  .  <N.>  +  <N.XN.>H.  .  (2-14) 


where  we  have  defined 


H.  .  -  -^  J  dld2  1^.(1,2)  (2-15) 

J    V  SI 

Now  because  of  the  translational  invariance  of  an  equilibrium 
ensemble  not  subject  to  external  fields 


hi.(l,2)  =  hij(R1,R2,  jflfn2)  =  hij(R1-R2,  n1ffi2)   (2"16) 

and  thus  we  can  also  write 

H.  .  =  -K-   J  dRdftndjLh.  .(0,,iL)  (2-17) 

13    vQ2    J      1   2  13  1'  1   2 

To  simplify  the  further  anaysis  we  rewrite  equation  2-14 
using  matrix  notation  as 

3<-N-  >  ,">     1  Q  \ 

f. i-)  =  [N  +  NHN].  .  (2-ia) 

where  the  elements  of  the  N  and  H  matrices  are 


N) .  .  =  6.  .  <N.>  (2-19a) 

'=  i]    °i]    l 
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<B>ij-8ij  (2"19b 


To  calculate  changes  in  chemical  potentials,  one  is 
interested  in  the  inverse  of  equation  2-18 


6^  r  "I, 

i_l  =  { [N  +  NHN]   } . .  (2-20) 


'3<V  T,V,<NV>         =    ===     ^ 

This  is  most  easily  expressed  in  terms  of  integrals 
of  the  direct  correlation  functions,  c.  .(1,2),  introduced 
by  Ornstein  and  Zernike  (1914).   These  direct  correlation 
functions  are  defined  by 


h.  j(l,2)  =  cij(l,2)  +  I  ^   J  d3cik(l,3)hkk(3,2)  (2-21 


Because  this  integral  is  not  of  full  convolution  form, 

it  may  seem  that  it  is  not  possible  to  relate  the  fluctuation 

derivatives  to  integrals  of  the  direct  correlation  functions 

but  that  is  incorrect.   Appendix  1  gives  the  details  of 

the  relationship  required.   Using  equations  2-20  and  Al-23 

we  find 

3Bu.  _1 

(^ M  =  [N    -  C/v]  .  .  (2-22) 

]   T,V,<N  > 

k  k^j 

and  if  the  independent  variables  used  are  the  component 
densities  the  result  is 
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Bu-  _i 

-1        =  [<L         ~    C].,  (2-23 


If  equation  2-23  is  combined  with  equation  2-10  one 
finds  that 

3  inf. 
f i]        =  -C.  .  (2-24) 

9pi   T  0  1D 

These  relations  are  most  useful  because  they  are 
required  to  derive  the  differential  equation  of  state. 
If  the  Gibbs-Duhem  equation  (at  constant  T)  is  written 


dPB  =  I    I    P±  [j^)  d  9.  <2"25) 

then  it  is  shown  that  the  equation  of  state  can  be  found 
through  knowledge  of  the  direct  correlation  function 
integrals.   Combining  equations  2-25  and  2-23  yields 

dPB  -  I  U  -  I  Pi  C^.ld  P,  (2-26) 

A  knowledge  of  the  C..  then  allows  for  the  calculation 
of  both  pressure  changes  and  fugacity  ratios  relative  to 
any  chosen  reference  state.   These  results  can  also  be 
used  to  determine  partial  molar  volumes  of  all  components 
in  a  mixture. 
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Summary 

This  chapter  has  dealt  with  some  of  the  properties 
of  interest  for  process  design.   The  relation,  required 
by  thermodynamic  consistency,  between  the  partial  molar 
volumes  and  the  chemical  potentials  has  been  emphasized. 
Mention  was  made  of  the  common  forms  of  correlations  for 
these  quantities,  asserting  that  often  liquid  reference 
state  approaches  are  used  for  both.   The  final  section 
presented  the  fluctuation  solution  theory  that  allows  for 
calculation  of  the  partial  molar  volumes  and  chemical  poten- 
tials, based  on  any  reference  state,  to  be  expressed  simply 
in  terms  of  ope  set  of  functions. 


CHAPTER  3 
GROUP  CONTRIBUTIONS 

One  of  the  more  powerful  tools  developed  for  physical 
property  estimation  has  been  the  group  contribution  concept. 
The  term  group  normally  refers  to  the  organic  and  inorganic 
radicals  but  can  be  more  specific.   A  molecule  of  interest 
can  be  described  by  the  number  of  the  different  types  of 
groups  of  which  it  is  composed. 

There  are  two  general  methods  for  using  the  group 
contribution  concept.   In  the  first  approach  some  molecular 
property  is  written  in  terms  of  the  state  variables  and 
a  set  of  parameters,  9, 

M  -  M(T,P;  9)  (3-1) 

and  the  parameters  for  a  given  set  of  substances  are  found 
as  sums  of  group  contribution 

9.  =    y    9.  (3-2) 

all 
groups  a 

The  most  common  forms  of  this  type  of  formulations  have 

been  for  predicting  critical  properties  (Lydersen,  1955) 

and  ideal  gas  specific  heats  (Verma  and  Doraiswamy,  1965). 
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In  some  cases  this  technique  can  be  justified  on  molecular, 
grounds . 

The  second  use  of  group  contribution  ideas  has  been 
to  assume  that  the  groups  actually  possess  thermodynamic 
properties  and  that  the  molecular  properties  are  then  a 
sum  of  these  group  properties, 

Mi  =     I        Ma  (3-3) 

all 
groups 
a  in  i 

This  idea  is  the  basis  for  two  popular  activity  coefficient 
correlations,  ASOG  (Derr  and  Deal,   1968)  and  UNIFAC 
(Fredendslund  et  al.,  1975).   A  model  for  the  group  proper- 
ties must  be  developed  for  this  technique  to  be  useful. 

The  true  utility  of  the  group  contribution  approach 
stems  from  its  predictive  ability.   Because  all  of  the 
molecules  in  a  homologous  series  are  formed  of  the  same 
groups,  only  in  different  proportions,  the  data  for  several 
of  the  elements  of  the  series  can  be  used  to  establish 
the  group  property  correlations.   These  can  then  be  used 
to  predict  the  properties  of  the  other  series  members. 
This  has  even  greater  scope  for  mixtures. 

Consider  for  example  mixtures  of  n-alkanes  and 
n-alkanols.   They  can  be  considered  to  be  made  of  only 
three  groups,  -OH,  -CH2,  and  -CH-j.   Thus,  any  mixture  of 
the  alkanes  and  alkanols  can  be  described  by  the  concentra- 
tion of  these  three  groups.   If  a  viable  theory  exists 
for  some  property  in  terms  of  the  group  functions,  then 
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the  properties  of  all  mixtures  of  these  groups  are  set. 
Figure  3-1  shows  an  application  of  these  ideas  for 
calculation  of  the  excess  enthalpy  of  alkane-alkanol  mix- 
tures.  The  figure  shows  the  surface  of  the  excess  enthalpy 
for  all  mixtures  of  the  hydroxyl,  methyl,  and  methylene 
groups  calculated  using  the  UNIFAC  equation  (Skjold- 
J^rgensen  et  al.,  1979).   Any  compound  made  of  the  three 
groups  is  represented  by  a  point  in  the  base  plane.   For 
example,  the  point  ( XCR   =  1/2,  XCH   =  0,  XQH  =  1/2)  is 
that  for  methanol.   The  possible  group  compositions  for 
any  mixture  are  found  along  the  line  connecting  two  molecular 
points.   In  the  figure  these  lines  are  drawn  in  for  methanol- 
pentane,  ethanol-pentane,  and  pentanol-pentane.   The  predic- 
tion of  the  excess  enthalpy  is  then  found  by  the  intersection 
of  a  vertical  plane  along  the  composition  path  and  the 
property  surface.   To  obtain  the  enthalpy  prediction  of 
the  molecular  system,  the  ideal  solution  value  must  be 
subtracted  from  the  group  estimate.   The  ideal  solution 
line  simply  connects  the  property  surface  at  the  points 
of  the  pair  molecules.   The  enthalpy  prediction  for  a 


ME 


molecular  system  of  methanol  and  pentane  and  at  (X 
2/3,  X    =  1/3)  is  shown  as  the  value  H   in  the  figure. 
The  power  here  is  that  this  one  diagram  can  be  used  to 
find  excess  enthalpies  for  all  alkane-alkanol  mixtures. 
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Fiaure  3-1.   Excess  enthalpy  surface  for  n-alkane-n-alkanol 
systems  at  298  K  calculated  using  UNIFAC 
theory. 
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"Reactive"  System  Theory 

Equation  3-3  can  also  be  written  as 

M   =  J"  v.  .   M  (3-4) 

i       in   a 

where  v.  ,  the  stoichiometric  coefficient,  represents  the 
number  of  groups  of  type  a  in  molecule  i.   This  expression 
is  completely  analogous  to  that  found  for  a  system  of  groups 
"reacting"  to  form  the  molecule 


7  v .  a 

L       ia 


(3-5 


With  this  idea  the  group  contribution  expressions  have 
a  physical  interpretation  and  the  thermodynamics  of 
"reactive"  systems  (Perry  et  al. ,  1981)  can  be  applied 
to  obtain  many  results. 

The  motivation  here  has  been  to  use  the  fluctuation 
solution  theory  in  terms  of  the  direct  correlation  function 
integrals.   It  is  then  desirable  to  determine  the  relation- 
ship between  the  molecular  and  group  direct  correlation 
functions.   This  has  been  accomplished  in  a  very  general 
fashion  by  Perry  (1980).   The  development  is  lengthy,  but 
the  salient  features  are  presented  below. 

The  most  important  aspect  of  this  approach  is  that 
even  if  the  groups  do  not  have  a  thermodynamics,  there 
are  still  well-defined  correlations  between  groups.   This 
allows  for  the  development  of  the  Kirkwood-Buf f  theory 
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in  terms  of  group  fluctuation,  with  the  constraints  offered 
by  equation  3-5.   The  result  is 

38P. 

j   T,p. 


10r)„         -  >sT  (a'"1  -£')  si^  (3-6, 


where  p'  is  the  matrix  of  group  densities  and  C'  is  a  matrix 
of  group  direct  correlation  function  integrals. 

While  this  theory  is  formally  exact  for  the  "reactive" 
system,  it  can  be  difficult  to  apply.   The  correlation 
function  integrals  contain  both  inter-  and  intramolecular 
contributions.   For  the  total  correlation  functions  these 
effects  can  be  separated,  but  no  known  analogous  result 
exists  for  the  direct  correlation  functions  ( Lowden  and 
Chandler,  1979).   This  is  a  problem  that  becomes  of  great 
importance  in  attempting  to  model  the  correlation  function 
integrals . 

RISM  Theory 

Chandler  and  Andersen  (1973)  have  formulated  a  molecular 
theory  for  hard  sphere  molecules  known  as  RISM.   The  mole- 
cules are  assumed  to  be  composed  of  overlapping  hard  spheres 
or  groups.   They  show  how  the  molecular  Ornstein-Zernike 
equation  2-18  can  be  reduced  to  a  group  form  with  explicit 
separation  of  the  intermolecular  and  intramolecular  correla- 
tions.  This  allows  them  to  define  group  direct  correlation 
function  integrals  that  have  only  intermolecular 
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contributions.   The  development  is  detailed  but  the  essential 
result  is 


-,(1,2)  -  II  vlfl  vj3  Ca6[r7,  x~)  (3-7) 


Note  here  that  the  molecular  correlation  function  has  orien- 
tation dependence  even  though  the  group  functions, 
cQg(r.,r.),  are  written  for  spherically  symmetric  inter- 
actions.  Chandler  and  Andersen  discuss  many  attributes 
of  these  group  functions  but  put  no  emphasis  on  the  thermo- 
dynamic ramifications  of  these  findings. 

Differences  between  RISM  Theory  and 
"Reactive"  Solution  Theory 

The  major  difference  between  the  RISM  theory  and  the 
"reactive"  solution  theory  is  the  nature  of  the  direct 
correlation  functions.   In  Perry's  formulation  we  have 


c   (r   )  =  direct  correlation  function        (3-1 
between  group  a  and  group  8 


and  the  RISM  theory  uses 


c   (ra,r  )  =  direct  correlation  function      (3-9) 
a®      i       J  between  group  a  on  molecule  i 

and  group  B  on  molecule  j 


This  shows  that  the  RISM  functions  contain  less  information 
but  may  be  easier  to  model  because  they  seem  analogous 
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to  molecular  functions  for  which  models  have  been  developed. 
Chandler  and  Andersen  have  also  presented  a  variational 
theory  which  can  be  used  to  obtain  the  direct  correlation 
functions  for  hard  body  systems. 

Figure  3-2  shows  the  results  obtained  using  data  from 
Lowden  and  Chandler  (1973)  for  a  system  of  hard  diatomic 
molecules.   The  diatomic  molecules  are  treated  as  overlapping 
spheres  of  diameter  o   separated  by  a  distance  L.   The  figure 
shows  that  the  calculations  agree  to  a  reasonable  degree 
with  the  Monte  Carlo  calculations  for  this  type  of  system. 

Another  important  aspect  of  the  Chandler-Andersen 
theory  is  that  they  discuss  how  the  cag  functions  could 
be  written  for  real  molecules.   In  general,  they  show  that 
perturbation  theories  that  would  be  valid  for  molecular 
systems  would  also  apply  to  RISM  group  system.   The  RISM 
formulation  will  be  employed  in  the  present  work. 

Thermodynamic  Properties  from  Group  Functions 

Equations  2-23  and  2-26  can  be  combined  with  equation 
3-7  to  express  the  molecular  thermodynamic  differentials 
in  terms  of  group  correlation  function  integrals.   The 
first  quantity  required  is 

C.  .  =  (-1?)    I  dld2  c.  .(1,2)  (3-10) 

Using  equation  3-7  this  is 
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Figure  3.2.  Comparison  of  RISM  ( )  theory  prediction  and  Monte 

Carlo  (•  •)  for  a  hard,  homonuclear  diatomic  molecule 
with  separation  to  diameter  ratio  of  0.6. 
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C-.  =  I  I    v.   v.Q  f-i-o)  I  dld2  c  R  (r\  rj)       (3-11 
i]    J     ia   D6  lv2fi2^  aB    i    D 


To  evaluate  this  quantity,  a  coordinate  transformation 
is  required. 

dld2  =  dR1dR2d^1,  d^2  ■*■   dr™,  dr^di^d^        (3-12) 

Here,  £2.,  represents  the  set  of  angles  needed  to  specify 
the  orientation  of  molecule  i  with  respect  to  a  fixed 
coordinate  system.   The  above  transformation  is  canonical, 
and  the  angular  integrations  can  be  performed  to  yield 

C..  =  Hv.   v.fl  [X)    J  d?ad?6  c  R(?a,  P.)  (3-13) 

a  B  V 

Now,  the  systems  under  consideration  are  homogeneous  so 
that 


c  R(?a,  ?S)  =  c  R(?a  -  ?6)  =  c   (?)  (3-14) 

aB   id      aB   i     3      aB 


Then 


C  •  =  I  I  v.   v..  C  Q  (3-15) 

13    L  „   ia   jB   aB 
a  B 


where 


C  Q  =  (i)  J  d?  c  _(?)  (3-16) 

aB    lV7  J      aB 


This  leads  to 
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=  r^2   ~    1    I    v.n    *,fl  CnR  (3-17a 


3p  •  ''  P  •         a        ia    J8   aB 


or  in  differential  form 


dp  . 

dBu.  =  — -  -  I    v.  ([C  .  dPft)  (3-17b) 

l    p  .  L      ia  q   aB    P 

l  a  B 


and  the  corresponding  expression  for  the  pressure  variation 


dPB  =  dP  -  I  I  Pa  C    dpB  (3-18) 

a  b 

Summary 

Group  contribution  approaches  are  valuable  for  predict- 
ing thermodynamic  properties.   In  this  chapter  two  general 
approaches  for  incorporating  group  contributions  into  fluc- 
tuation theory  have  been  presented  and  analyzed.   Finally, 
the  thermodynamic  property  differentials  for  molecular 
systems  have  been  expressed  in  terms  of  group  direct  correla- 
tion function  integrals.   Actual  property  changes  can  be 
formulated  once  models  are  expressed  for  the  integrals. 
This  is  the  subject  of  the  next  chapter. 


CHAPER  4 
FORMULATION  OF  MODELS  FOR 
DIRECT  CORRELATION  FUNCTIONS 


In  previous  chapters  it  has  been  shown  that  changes 
in  thermodynamic  properties  of  molecular  systems  can  be 
expressed  in  terms  of  integrals  of  site  direct  correlation 
function.   Expressions  for  these  integrals  are  required 
before  the  calculations  can  be  performed. 

The  purpose  of  this  chapter  is  two-fold.   First,  thermo- 
dynamic requirements  on  models  for  the  direct  correlation 
functions  will  be  presented  and  examined.   With  these 
restrictions  on  model  form  established,  a  feasible  model 
for  the  correlation  function  integrals  will  be  presented. 

Thermodynamic  Consistency 

One  aspect  that  has  been  emphasized  throughout  this 
work  is  the  necessity  of  the  formalism  to  meet  thermodynamic 
consistency  requirements.   Because  our  proposed  calculational 
scheme  employs  direct  correlation  function  integrals,  the 
models  that  are  developed  for  these  quantities  are  subject 
to  consistency  tests.   The  easiest  check  that  can  be  employed 
is  one  of  equality  of  cross  partial  derivatives  of  the 
(dimensionless )  chemical  potentials. 
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)p  .  dp,  • 

]   k 
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3p,  3p  . 
k   j 


V  i,  j,k 


(4-1) 


This  can  be  expressed  in  terms  of  the  direct  correlation 
functions  integrals  because 


36yi 


i:  _ 


3p^   *-pWj     0i 


ID 


(4-2) 


Combination  of  equations  4-1  and  4-2  leads  to 


3C 


3C 


ik- 


k   T,p 


l^k 


j   T,p 


(4-3) 


£^j 


It  must  be  noted  that  this  is  equivalent  to  the  equality 

of  three-body  direct  correlation  function  integrals  (Brelvi, 

1973) 


C  .,  =  C,  . 
ljk    ik] 


(4-4 


where  the  subscripts  can  take  on  all  values  associated 
with  the  species  in  the  mixture.   The  case  of  interest 
here  is  that  in  which  the  direct  correlation  function 
integrals  for  the  molecular  species  are  written  in  terms 
of  group  quantities 


Cij  =  l\    Via  Vj6  CaB 

a  S 


Then  equation  4-3  takes  the  form 


4-5 
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^7  [l    I    Via  Vj3  Ca6]^ 
k   a  6 


=   Clf:  II  J  yia  vkB  CaB  5  ^-6) 

This  requirement  can  be  expressed  in  terms  of  group  prop- 
erties if  the  chain  rule  is  used  for  the  derivatives 


af"  =  I     ^  W~  =  l   viy  if"  (4"7) 

i    y      i     Y    Y        Y 


Combination  of  equations  4-6  and  4-7  leads  to 
III    "*»  >V^  '^'t.P^ 

"     III  W*^,.,,^         «- 

This  relation  will  only  be  satisifed  in  general  if 

3C  .  3C 

( ajn        _  i ayi 


,pY   T  o        9pB   T  d  C4_9) 

Y   T,pn^Y       B   T,pn^6 


Thus,  the  consistency  requirements  for  group  direct  correla- 
tion function  integrals  are  equivalent  to  those  of  the 
molecular  quantities.   This  constraint  will  be  employed 
in  formulation  of  the  working  models  for  the  group  direct 
correlation  function  integrals. 
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Basis  for  Model  Development 

It  has  already  been  shown  that  the  thermodynamic  con- 
sistency tests  for  a  group  formulation  and  a  molecular 
formulation  are  equivalent  when  written  in  terms  of  direct 
correlation  function  integrals  based  on  the  RISM  theory. 
This  should  not  be  surprising  considering  that  stability 
conditions  for  fluids  were  shown  to  be  equivalent  for  the 
two  approaches  by  Perry  (1980).   Thus,  in  developing  models 
for  group  direct  correlation  function  integrals,  the  same 
considerations  should  apply  as  those  used  by  Mathias  (1979) 
in  developing  molecular  models. 

The  general  philosophy  that  will  be  employed  is  to 
use  as  much  theoretical  information  as  possible  to  develop 
these  models.   This  requires  use  of  some  concepts  of  statis- 
tical mechanical  perturbation  theory  to  obtain  approximate 
forms  for  the  correlation  function  integrals.   The  analysis 
of  Chandler  and  Andersen  (1972)  has  shown  that  for  RISM 
theory  direct  correlation  function  results  of  molecular 
perturbation  theory  are  easily  extended  to  the  group  func- 
tions.  Appendix  2  contains  a  complete  development  of  an 
exact  perturbation  theory  for  the  direct  correlation  func- 
tions based  on  the  RISM  theory.   In  this  chapter  we  shall 
only  deal  with  the  important  ramifications  of  these  results. 
It  is  always  possible  to  write 


_  rref    +  (r        -  cref )  (4-10) 

-a  6  "  CaS   +  (Ca8    Ca  6  '  ^ 
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where  the  superscript  ref  refers  to  some  reference  system. 
The  purpose  of  perturbation  theory  is  to  determine  a  refer- 
ence system  so  that  an  approximate  form  for  the  perturbation 
(second  on  the  right-hand  side)  term  can  be  made  that  yields 
useful  results.   Because  this  work  is  concerned  with  dense 
fluids,  we  require  a  reference  system  that  can  adequately 
represent  the  behavior  at'  high  densities.   In  liquids  the 
optimal  choice  of  a  reference  system  seems  to  be  one  with 
only  repulsive  forces  (Weeks,  Chandler,  and  Andersen,  1971). 
Values  of  C  g  are  not  available  for  this  type  of  model 
system.   However,  the  system  with  purely  repulsive  forces 
can  be  well  represented  by  a  system  of  hard  spheres  if 
the  hard  sphere  diameters  are  chosen  as  functions  of  tempera- 
ture (Barker  and  Henderson,  1967).   This  approach  will 
be  followed  in  this  work.   Even  with  this  choice  of  reference 
system  the  perturbation  term  cannot  be  exactly  identified. 
A  further  approximation  used  is  that  the  zero  density  limit 
of  this  term  is  adequate.   Thus,  the  model  employed  here 


r         -    rHS    +   lim  (c         -    CEh  (4-11 

caS  ~  Ca6  +   p+o  (Ca8    La6J  { 


Appendix  2  shows  the  evaluation  of  the  required  limit  which 
gives  the  working  relation 


C  ,  =  CH?  +   I   (a  )   /Tn  (4-12) 

aB     aB    n=o    n  aB 
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The  constants  in  equation  4-12  are  dependent  on  the  inter- 
molecular  potential  (written  as  a  sum  of  group  of  interac- 
tions) .   No  explicit  calculation  of  these  terms  is  attempted 
here  because  we  choose  to  determine  the  expressions  on 
the  basis  of  experimental  data.   This  model  form  is  com- 
pletely analogous  to  the  model  used  by  Mathias  (1979). 

Analysis  of  the  Model 

Several  interesting  features  of  the  direct  correlation 
function  model  seem  to  merit  mention.   Equation  4-12  is 
highly  similar  to  the  RISM  form  of  the  mean  spherical 
approximation  (Lebowitz  and  Percus,  1966;  Chandler  and 
Andersen,  1972)  and  is  essentially  equivalent  if  the  hard 
sphere  diameters  are  chosen  as  functions  of  temperature 
only.   The  assumptions  involved  in  obtaining  the  present 
approximant  and  the  mean  spherical  form  are  different, 
but  it  seems  that  for  our  purposes  this  difference  is 
immaterial. 

Of  greater  interest  here  is  the  relationship  suggested 
by  the  present  model  for  the  three-body  direct  correlation 
function  integrals.   These  can  be  found  from 

3C 
C.  ..  =  i^r2)  (4-13) 

The  form  for  the  C.  .  proposed  in  equation  4-12  then  suggests 

that 
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C.  ..   *  CHSV  (4-14) 

ljk    ljk 


This  is  probably  not  a  very  accurate  approximation,  but 
the  flexibility  in  the  model  form  inherent  due  to  the 
determination  of  the  parameters  from  experimental  data 
may  make  this  adequate.   A  better  interpretation  of  this 
analysis  is  that  the  density  dependence  of  the  three-body 
direct  correlation  function  is  assumed  to  be  approximated 
by  that  of  the  hard  sphere  quantity.   This  analysis  also 
suggests  an  extension  of  the  RISM  theory  to  three-body 
correlation  functions,  a  derivation  of  which  can  be  found 
in  Appendix  3. 


C     =JJJv.v.0vCo  ( 4-15 

ijk   £  §  ^   i«  jS  ky  S*By 


Finally,  thermodynamic  consistency  of  the  proposed  model 
must  be  considered.   If  equation  4-12  is  rewritten  in 
explicit  form 


CaS  <£rT>  =  Ca6  (^'T)  +  Kr    (T)  (4"16 


then  the  consistency  requirement  is  met  if 
SC^  3CHS 

3PY   T,p   ,       3pB   T,p   ,. 


4-17 


In  this  work,  the  C^p  are  determined  from  the  generalized 
hard  sphere  equation  of  state,  which  is  consistent,  then 
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this  must  be  true.   Note  that  if  terms  of  higher  order 
in  density  had  been  included  in  the  perturbation  term, 
they  also  would  have  to  have  been  chosen  to  be  thermodynami- 
cally  consistent. 

Also  of  interest  here  is  the  functional  dependence 
of  the  hard  sphere  diameters.   If  they  are  chosen  to  be 
functions  of  density,  as  is  normally  done  in  the  Week- 
Chandler-Andersen  perturbation  scheme,  then  this  dependence 
would  have  to  be  such  that  equation  4-17  was  satisfied. 
This  is  an  aspect  not  always  appreciated. 

Property  Changes  Using  the  Proposed  Model 

This  chapter  has  presented  and  analyzed  a  feasible 
model  for  group  direct  correlation  function  integrals. 
The  utility  of  this  model  can  only  be  tested  by  its  use 
to  calculate  changes  in  molecular  thermodynamic  properties. 
The  required  formulas  are  developed  below. 

Chemical  Potential 

Equation  3-17b  expresses  the  differential  of  a  molecular 
chemical  potential  in  terms  of  the  group  direct  correlation 
function  integrals.   This  is  conveniently  rewritten  in 
terms  of  the  residual  chemical  potential  as 


dS^  .  I    I    V    caf)  dP6  (4-18) 

a  p 


Insertion  of  equation  4-16  and  integration  leads  to 
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«Mj    -J    vla     [ABPar'HS    +    [    *„B    4p6I  (4-19) 

a  p 

where  ^£'HS  is  the  residual  chemical  potential  calculated 
from  the  hard  sphere  equations  for  group  a  in  a  solution 
of  groups. 

Pressure  Change 

To  calculate  pressure  changes,  it  is  again  easiest 
to  work  with  residual  properties.   The  required  differential 
f o rm  is 


dPDr  =  "  I  I  P  C  R  dP 
8       a  3 


4-20 


Insertion  of  the  model  and  integration  leads  to 

APB   =  APB'HS  +  <!>  H    Kt      ^Pap8)  (4_21) 

where  Ap„'HS  is  the  change  in  the  residual  pressure  (divided 

p 

by  k  T)  calculated  from  the  hard  sphere  equation  for  the 
B 

solution  of  groups. 

It  should  be  noted  that  the  ideal  gas  state  used  for 
the  residual  property  changes  on  either  side  of  equation 
4-21  are  different.   This  should  be  expected  for  the  group 
equation  of  state  is  written  in  terms  of  the  group  densities 
which  sum  to  a  larger  value  than  the  molecular  densities. 
However,  beacause  this  form  is  not  used  as  a  complete  equa- 
tion of  state,  this  should  not  present  a  problem. 


CHAPTER  5 
MODEL  PARAMETERIZATION 


This  work  has  been  concerned  with  formulating  the 
thermodynamic  properties  of  molecular  systems  in  terms 
of  integrals  of  group  direct  correlation  functions.   In 
the  last  chapter  a  form  for  these  terms  was  developed 


_  rREF    ^  (5_1 


and  further,  the  reference  state  was  identified  as  that 
of  a  hard  sphere  system.   The  full  dependence  on  state 
variables  is  written  as 


Ca6  (^'T)  =  CS(3  (-'T;  -]    +    ^aB  (1/T)  (5~2 


where  _£  is  the  vector  of  hard  sphere  diameters  for  the 
groups  in  the  system.   The  purpose  of  this  chapter  is  to 
develop  empirical  forms  for  the  functions  ^a„  and  to  estab- 
lish the  dependence  of  the  hard  sphere  diameters  on  state 
parameters  (Weeks,  Chandler,  and  Andersen,  1971). 
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Corresponding  States  Theory 

The  utility  of  the  form  proposed  in  equation  5-2  can 
be  enhanced  if  the  perturbation  functions  and  hard  sphere 
diameters  could  be  written  in  corresponding  states  form. 
This  is  a  two-step  process: 

1.  Determine  if  a  corresponding  states  principle 
exists  for  the  direct  correlation  function 
integrals  for  liquids. 

2.  Determine  if  the  functional  form  proposed  in 
equation  5-2  can  represent  the  observed 
correspondence. 

To  answer  these  questions  simple  molecules  will  be 
considered  first.   The  molecules  considered  here  can  be 
considered  to  consist  of  only  one  group.   The  relation 
to  thermodynamics  used  in  the  analysis  is 


Hj)   .  1/U-PC) 


Figure  5-1  shows  a  corresponding  states  correlation  for 

the  quantity  on  the  left-hand  side  of  equation  5-3  for 

simple  molecules,  using  only  one  parameter.   Further,  Brelvi 

and  O'Connell  (1972)  showed  that  this  behavior  can  be  found 

for  polyatomic  molecules  also,  if  the  density  is  large 

enough  (  P  >  2  p  ) .   These  results  suggest  that  for  large 
r  c 

densities  a  one-parameter  corresponding  states  formulation 
may  be  valid  for  the  direct  correlation  function  integrals. 
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Figure  5-1.   Bulk  modulus  versus  reduced  ensity  for  simple 
fluids.   The   line  is  argon  data  and  0  is 
Lennard-Jones  potential  parameter. 
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Note  that  the  temperature  dependence  seems  to  have  little 
effect.   Inclusion  of  temperature  dependence  in  these  cor- 
relations should  be  able  to  enhance  the  accuracy. 

Figure  5-2  shows  another  test  of  this  correspondence. 
This  is  in  an  integrated  form  of  equation  5-3. 

-     P 
P6  =  P8rer  +  J  (1-PC)  dP  (5-4) 

pref 

Here,  a  characteristic  volume  for  each  molecule,  V*,  and 
a  characteristic  temperature,  T*,  have  been  found  to  yield 
the  correspondence  as 


P  =  PrSf'+  f(p,  T;  p~ref)  (5-5) 


where 

P  =  PV*/RT 

p  =  pV* 

T  =  T/T* 

The  function,  f,  in  equation  5-5  is  a  corresponding 
states  correlation  for  the  density  integral  of  the  direct 
correlation  function.   This  is  what  is  desired.   The  data 
presented  seem  to  suggest  that  a  corresponding  states  formu- 
lation can  be  found  for  pure  component  direct  correlation 
function  integrals  for  dense  fluids.   A  detailed  discussion 
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Figure  5-2.   Corresponding  states  correlation  for  liquids 
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of  the  origins  of  this  correlation  is  presented  in  Appendix 

4. 

Mathias  (1979)  has  shown  that  these  concepts  can  be 
extended  to  mixture-direct  correlation  function  integrals 
with  the  proper  choice  of  mixing  rules.   He  used  a  form 
similar  to  equation  5-2,  but  the  group  contribution  concept 
was  not  employed. 

The  applicability  of  the  two-parameter  corresponding 
states  principle  seems  valid  for  the  direct  correlation 
function  integrals.   The  question  of  the  use  of  the  func- 
tional form  must  now  be  addressed.   This  is  again  most 
easily  accomplished  by  analysis  of  a  set  of  experimental 
data  for  a  simple  molecule.   After  the  noble  gases  are 
used  for  this  type  of  testing,  however,  the  major  interest 
here  is  for  hydrocarbons  so  that  a  hydrocarbon  reference 
system  may  prove  to  be  more  useful. 

Choice  of  a  Reference  Component 

Mathias  (1979)  had  great  success  in  developing  a 
corresponding  states  correlation  for  direct  correlation 
function  integrals,  using  a  form  similar  to  equation  5-2, 
for  molecular  systems,  using  argon  as  a  reference  component 
to  determine  model  parameters.   This  work  is  concerned 
with  hydrocarbon  properties  so  that  it  seems  reasonable 
to  examine  the  difference  in  direct  correlation  function 
integrals  between  argon  and  a  simple  hydrocarbon,  methane. 
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The  comparison  is  performed  here  using  the  dimensionless 
quantity  1-pC.   Table  5-1  shows  this  comparison  for  two 
isotherms,  in  the  dense  fluid  region. 

The  behavior  is  similar  on  both  isotherms.   It  is 
seen  that  the  fluids  act  quite  similarly  at  lower  densities 
but  behave  differently  as  the  density  becomes  larger. 
Methane  is  increasingly  more  compressible  as  the  density 
increases.   This  can  be  explained  by  the  polyatomicity 
of  methane.   Both  molecules  are  essentially  spherical, 
but  at  the  higher  densities  the  methane  seems  to  allow 
for  interlocking  of  the  hydrogens. 

These  effects  can  be  important  for  the  development 
of  the  group  contribution  correlations.   If  the  groups 
were  chosen  as  monatomic,  then  argon  would  probably  be 
the  better  choice  of  reference  substance  due  to  physical 
similarity.   However,  if  the  groups  are  chosen  as  the  common 
organic  radicals,  then  methane  may  be  the  better  choice 
of  reference  substance. 

The  actual  choice  of  a  reference  substance  will  be 
made  after  it  has  been  determined  whether  equation  5-2 
is  a  viable  functional  form  for  the  direct  correlation 
function  integrals. 

Method  of  Data  Analysis 

This  section  addresses  the  question  of  whether  the 
form  proposed  in  equation  5-2  can  be  used  to  fit  the  direct 
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TABLE  5-1 

COMPARISON  OF  DIRECT  CORRELATION  FUNCTION  INTEGRALS 
FOR  ARGON  AND  METHANE 


T/Tc  = 

0.76 

T/T   = 
/  c 

2.0 

pvc 

(1-PC) 

argon 

(  1- 

p   methane 

pvc 

(1-pC) 

argon 

(1-f 

3   methane 

2.26 

6.44 

6.35 

2.1 

7.23 

7.14 

2.31 

7.97 

7.82 

2.2 

8.50 

8.35 

2.37 

9.68 

9.44 

2.3 

9.96 

9.73 

2.43 

11.51 

11.25 

2.4 

11.  65 

11.29 

2.48 

13.57 

13.24 

2.5 

13.58 

13.07 

2.54 

15.83 

15.43 

2.6 

15.79 

15.06 

2.59 

18.30 

17.83 

2.7 

18.29 

17.30 

2.65 

21.05 

20.46 

2.8 

21.  13 

19.91 

2.70 

23.88 

23.33 

2.9 

24.33 

22.59 

2.76 

27.10 

26.45 

3.0 

27.94 

25.68 
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correlation  function  integrals  of  real  fluids.   The  analysis 
will  be  performed  in  terms  of  the  dimensionless  quantity, 
C ,  defined,  by 

3PB 
C  S  PC  -  1  -  (v-)t  (5-6. 

The  thermodynamic  derivative  is  evaluated  from  correlated 
experimental  data.   From  equation  (5-2)  this  is 

C  =  CHS  (p,  T;  a)    +  pip  (5-7) 

For  this  work  the  hard  sphere  diameter  is  treated  as  a 
function  only  of  temperature.   Thus,  along  an  isotherm, 
if  the  proposed  functional  form  is  proper,  for  the  proper 
choice  of  hard  sphere  diameter,  the  value  of  i>    should  be 
constant  given  by 

<|>  =  (cHS  -  C)/p  (5-8) 

for  any  value  of  p.   The  suitability  of  the  proposed  model 
is  determined  by  how  well  this  is  obeyed.   The  hard  sphere 
properties  are  determined  using  the  generalized  hard  sphere 
expressions  as  presented  in  Appendix  5.   It  should  be  noted 
that  all  of  the  hard  sphere  equations  are  generated  by 
the  same  microscopic  form  of  the  direct  correlation 
function.   Here  the  hard  sphere  direct  correlation  function 
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integrals  are  found  using  equation  5-6  with  the  corresponding 
form  of  the  hard  sphere  equation  of  state.   For  these  pure 
component  cases  we  have 


4  3      2 

4(1+1     '   1" 

(1-n) 


-hs  _  n+a)n  -  4(l+a)n  +  2i ~_ 

C    —  4 


3 
where   n  =  n/6.  po 


and  in  corresponding  states  form 

n=    [  tt/6    o3/V*][pV*]  (5-10a) 

e   X    pV*    =    Xp  (5-10b) 

Here,  the  characteristic  volume,  V*,  has  been  introduced 
as  a  reducing  parameter  for  both  the  hard  sphere  volume 
and  the  density.   On  each  isotherm  there  are  actually  two 
parameters  that  can  be  used  to  vary  the  model,  a  and  X. 
The  analysis  here  has  been  done  sequentially,  a  value  of 
a  is  chosen  and  then  the  optimum  X  has  been  found  for  that 
value  of  a.   The  characteristic  volume  has  also  been  used 
to  reduce  the  perturbation  term  i|>  as 

p  =  (ii>/v*)V*  =  V*f(T)  (5-11) 

The  model  can  then  be  written  as 
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C(P,T)  =  CHS(n;  T,X,a)  +  pf  (5-12) 

Figure  5-3  gives  a  simplified  flowchart  of  how  the  calcula- 
tions are  performed  for  each  isotherm  for  a  fixed  value 
of  a.   The  subscript  i  denotes  the  number  of  data  points 
on  the  isotherm,  NP.   The  optimization  routine,  RQUADD, 
can  be  found  in  Appendix  6. 

Analysis  of  Argon  Data 

3PB 
The  values  of  (-= — L  needed  for  this  analysis  were  generated 

v  dp   'T 

from  the  equation  of  state  of  Twu   et  al.  (1980).   These 
authors  claim  to  be  able  to  reproduce  dense  fluid  pressures 
to  within  the  experimental  accuracy.   Eleven  temperatures, 
evenly  spaced  between  the  triple  and  critical  points,  were 
chosen  to  generate  liquid  phase  data.   Each  isotherm  con- 
sisted of  16  density-correlation  function  integral  pairs 
evenly  spaced  between  the  vaporizing  and  freezing  densities. 
Supercritical  data  were  also  used.   The  isotherms  were 
for  values  of  T/T   from  1.1  to  3.0.   The  densities  used 
were  in  0 . 1  increments  of  P /P   from  1.5  to  3.0,  or  the 
freezing  density,  whichever  was  lower.   A  listing  of  the 
program  used  to  generate  the  data  is  given  in  Appendix 
6. 
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Input  X 


find  f.  =  [C    (P.)  -  CHC(P.  ,ot)  ]/P. 
1      exp   1      HS   l       l 


NP 
find  f  =  ^p  I      fL 


update  X  using ^ RQUADD 


find  C^alc  =  CH(Pif  a)  +  p±  f 


NP 

n  «,^     V  I  ^calc    ~exP 

find  AAD  =   L  I  C.     -  C.  ^ 

i=l    x 


is  AAD.  a  minimum 


print  results 


Fiqure  5-3.   Flowchart  for  calculation  of  X  at  a  fixed  a. 
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For  the  initial  analysis  values  of  a  of  0,  -1,  -2, 
-3,  -4,  and  -5  were  used.   Table  5-2  shows  the  average 
absolute  deviation  in  the  direct  correlation  function 
integrals  found  for  the  optimum  value  of  x  at  each  isotherm. 
A  noticeable  minimum  exists  in  the  range  of  a  =  -4.   Note 
that  the  Carnoban-Starling  equation  (a  =  -1)  does  not  yield 
the  best  predictions,  even  though  it  is  the  "best"  hard 
sphere  equation  of  state. 

Several  isotherms  were  chosen  to  examine  the  sensitivity 
of  the  results  to  the  value  of  a,  on  a  finer  scale.   These 
results  are  given  in  Table  5-3.   For  the  liquid  phase  iso- 
therms a  value  of  a  =  -4.3  seems  to  be  optimal.   The  super- 
critical data  are  best  fit  with  a  of  -4.0  or  -4.1.   It 
was  decided  to  use  a  constant  value  of  a  =  -4.2  to  fit 
all  of  the  argon  data. 

With  a  set  at  -4.2  the  value  of  X  and  f  could  then 
be  determined.   The  optimum  values  of  X  found  were  correlated 
as  a  polynomial  in  inverse  temperature.   For  this  analysis 
the  characteristic  parameters  were  chosen  to  be  equal  to 
the  critical  parameters.   The  form  that  fit  X  best  was 


0  12342  +  0-15069  _  0-18545  +  0.15947 

X  *  -v  _ 

i  T  T 


0.070261  0.012111  (r     nl 

^4  +  ^5  [  ' 

T  T 
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TABLE  5-2 
REGRESSION  ANALYSIS  OF  ARGON  DATA 


T 

Average  Absolute  Deviation  in  C  for 

=  a 

0 

-1 

-3 

-4 

-5 

0.5966 

0. 

0101 

0.0102 

0.0094 

0.0058 

0. 

7266 

0.6297 

0. 

0392 

0.0395 

0.0361 

0.0219 

1. 

3274 

0.6629 

0. 

0785 

0.0789 

0.0714 

0.0427 

1. 

7500 

0.6960 

0. 

1225 

0.1228 

0.1100 

0.0642 

2. 

0472 

0.7292 

0. 

1651 

0.1651 

0.1461 

0.0830 

2. 

2372 

0.7623 

0. 

2062 

0.2059 

0.1798 

0.0983 



0.7954 

0. 

2445 

0.2436 

0.2094 

0.1085 

2. 

3746 

0.8286 

0. 

2758 

0.2736 

0.2304 

0.1104 

2 

3629 

0.8949 

0. 

3104 

0.3056 

0.2408 

0.1122 

2 

1992 

0.9280 

0 

3053 

0.2990 

0.2374 

0.1280 

2 

0497 

0.8617 

0 

2994 

0.2963 

0.2430 

0.1071 

2 

2885 

1.1 

0 

2840 

0.2764 

0.2142 



1 

4495 

1.2 

0 

2413 

0.2344 

0.1835 

0.1260 

1 

1049 

1.3 

0 

2175 

0.2088 

0.1580 

0.1058 

0 

8986 

1.4 

0 

1990 

0.1899 

0.1383 

0.0912 



1.5 

0 

1837 

0.1740 

0.1227 

0.0801 

0 

.6364 

1.6 

0 

.1707 

0.1611 

0.1101 

0.0713 

0 

.5516 

1.7 

0 

.1611 

0.1504 

0.1003 

0.0643 

0 

.4845 

1.8 

0 

.1530 

0.1411 

0.0920 





1.  9 

0 

.1498 

0.1344 

0.0854 

0.0536 

0 

.3864 

2.0 

0 

.1403 

0.1286 

0.0798 

0.0495 

0 

.3505 

2.2 

0 

.1310 

0.1196 

0.0713 

0.0425 

0 

.2948 

2.4 

0 

.1242 

0.1131 

0.0651 

0.0371 

0 

.2536 

2.6 

0 

.1192 

0.1084 

0.0610 

0.0324 

0 

.2214 

2.8 

0 

.1152 

0.1047 

0.0582 

0.0290 

0 

.1960 

3.0 

0 

.1126 

0.1023 

0.0566 

0.0261 

0 

.1749 
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TABLE  5-3 
EFFECT  OF  HS  EQUATION  ON  FITTING  ARGON  DATA 


T 

Avera 

ge  Abso 

lute  Deviation  in 

C  for 

=  a 

-3.8 

-3.9 

-4.0 

-4.1 

-4.2 

-4.2 

-4.4 

.7954 

.1435 

.1277 

.1085 

.0848 

.0547 

.0214 

.0401 

.8286 

.1516 

.1329 

.1104 

.0831 

.0535 

.0263 

.0595 

.8617 

.1507 

.1297 

.1071 

.0829 

.0592 

.0399 

.0866 

.8949 

.1477 

.1294 

.1122 

.0916 

.0732 

.0667 

.1264 

.9280 

.1581 

.1434 

.1280 

.1144 

.1049 

.1152 

.1851 

1.1 

.1484 

.1368 

.1309 

.1235 

.1258 

1.2 

.1344 

.1303 

.1260 

.1263 

.1364 

1.3 

.1126 

.1092 

.1058 

.1064 

.1065 

1.8 

.0607 

.0590 

.0585 

.0618 

.0698 
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The  maximum  error  in  calculated  X  values  from  equation 
5-13  was  0.105%  with  an  average  error  0.022%.   With  X  values 
from  equation  5-13  and  a  set  at  -4.2  the  optimal  values 
of  f  were  then  found.   These  values  were  also  fit  to  a 
polynomial  as 


14.70    23.237  _,_  25.221    11.636 

f  =  -1.4098  + Z2 +  ^3 ^4 

T        T         T        T 


2.0463  (5_14 

T5 


This  polynomial  form  reproduced  the  f  values  with  an  average 
error  of  0.15%  for  the  isotherms  analyzed.   The  temperature 
dependence  of  X  and  f  is  shown  in  Figure  5-4.   Table  5-4 
gives  a  summary  of  the  values  of  X  and  f  for  argon  along 
with  the  average  error  in  the  prediction  of  the  direct 
correlation  function  integrals.   In  general,  the  correlation 
function  integrals  were  correlated  to  within  the  experimental 
accuracy.   The  only  range  of  temperatures  for  which  the 
fit  was  not  excellent  was  near  the  critical  point,  which 
is  to  be  expected. 

Analysis  of  Methane  Data 

There  is  less  high  quality  P-V-T  data  available  for 
dense  methane  than  for  argon.   In  this  work  the  equation 
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Figure  5-4.   Dimensionless  hard  sphere  diameter  and 
perturbation  function  based  on  argon 
reference . 


TABLE  5-4 
CORRELATION  OF  ARGON  DCFI  USING  a    =  -4.2 
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Average  Error 

T 

X 

f  V 
c 

in  C 

0.5966 

0.21159 

890.79 

.0035 

0.6297 

0.20916 

820.09 

.0128 

0.6629 

0.20692 

759.79 

.0525 

0.6960 

0.20482 

707.47 

.0372 

0.7292 

0.20283 

661.13 

.0466 

0.7623 

0.20094 

620.05 

.0532 

0.7954 

0.19915 

583.39 

.0562 

0.8286 

0.19745 

550.40 

.0547 

0.8617 

0.19583 

520.61 

.0596 

0.8949 

0.19430 

493.31 

.0732 

0.9280 

0.19286 

467.98 

.1076 

1.1 

0.18649 

373.74 

.1264 

1.2 

0.18348 

335.43 

.1377 

1.3 

0.18086 

306.28 

.1166 

1.4 

0.17852 

282.27 

.1037 

1.5 

0.17643 

262.19 

.0936 

1.6 

0.17453 

245.04 

.0844 

1.7 

0.17280 

230.17 

.0764 

1.8 

0.17120 

217.10 

.0700 

1.9 

0.16971 

205.48 

.0648 

2.0 

0.16833 

195.05 

.0605 

2.2 

0.16582 

177.04 

.0527 

2.4 

0.16359 

161.94 

.0458 

2.6 

0.16159 

149.05 

.0393 

2.8 

0.15978 

137.88 

.0343 

3.0 

0.15814 

128.10 

.0308 
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of  state  of  Mollerup  (1980)  was  used  to  generate  nine  sub- 
critical  isotherms  and  one  supercritical  isotherm,  as  was 
done  for  argon.   These  supercritical  isotherms  measured 
by  Robertson  and  Babb  (1969)  were  also  analyzed.   The  pro- 
cedure used  for  the  analysis  of  the  methane  data  was  the 
same  as  for  argon,  except  some  positive  values  of  a  were 
also  investigated. 

For  the  liquid  phase  isotherms  a  value  of  a  =  -4.2 
was  clearly  the  best.   For  the  supercritical  data  values 
of  a  from  -3.0  to  -4.3  could  produce  almost  equivalent 
results,  but  a  =  -4.0  was  the  optimum.   It  should  be  noted 
that  for  the  supercritical  isotherm  values  of  a  =  9  could 
also  offer  a  reasonable  representation  of  the  data,  but 
not  as  good  as  -4.0.   Thus,  a  value  of  -4.2  was  chosen 
to  represent  the  methane  data. 

Proceeding  as  was  done  before  for  the  argon  analysis, 
the  values  of  X  and  f  were  fit  to  polynomials  in  inverse 
temperature.   The  resultant  forms  were 


0.026446    0.080054    0.089994 

X  =  0.13832  +  z +  ^9 ^3 

T  T  T 


0.039085  _  0.0060984  (5-15a) 

T  T 


23.21         50.827         40.662         15.811 

f    =    5.0851 = +   r? =1 +   ^4— 

T  t  T  T 


2.3173 

-5 

T 


(5-15b) 
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where  the  critical  properties  were  used  as  the  reducing 
parameters.   These  functions  are  plotted  in  Figure  5-5. 
Table  5-5  gives  the  summary  of  results  for  the  methane 
correlation.   The  proposed  model  was  able  to  reproduce 
the  data  to  within  the  exerimental  uncertainty. 

Use  of  the  Correlations  to  Calculate  Pressure  Changes 

The  development  of  the  DCFI  models  for  the  pure  com- 
ponents was  first  tested  by  performing  equation  of  state 
calculations  for  the  base  substances,  argon  and  methane. 
The  method  of  calculation  was  to  choose  the  lowest  density 
point  on  each  isotherm  as  a  reference  and  then  to  calculate 
pressure  changes  for  all  other  points  on  the  isotherm  from 
that  reference.   There  are  thus  two  criteria  for  goodness 
of  fit,  the  average  percent  error  in  pressure  changes, 
defined  by 

ND     APCalC  -   Apf  ? 

AP  =  V  [  I  — i —  I  ]  /  ND  (5-16) 

H   '       Apf P 

and  the  sum  squared  percent  error,  defined  by 

ND   APCalc  -  APexP   2 
SSEAP  =  I    [         X   exp ±—    ]  (5-17) 

i=l      APSXP 

l 

where  in  both  cases  ND  is  the  total  number  of  data  points 
for  a  compound.  When  parameters  were  found  to  best  "fit" 
a  set  of  data,  the  quantity  SSEAP  was  minimized  using  a 
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TABLE  5-5 
CORRELATION  OF  METHANE  DCFI  USING  a  =  -4.2 


Average  Error 

T 

X 

f  V 
c 

in  C 

.4986 

.22182  . 

1545.26 

.0015 

.5511 

.21596 

1309.43 

.0155 

.6036 

.21095 

1131.07 

.0188 

.6561 

.20671 

996.35 

.0157 

.7085 

.20302 

890.93 

.0140 

.7610 

.19974 

805.32 

.0335 

.8135 

.19675 

733.53 

.0612 

.8660 

.19400 

671.78 

.0946 

.9185 

.19145 

617.18 

.1413 

1.1 

.18381 

473.60 

.1484 

1.6173 

.16916 

260.80 

.2262 

1.9585 

.16316 

204.86 

.1437 

2.4833 

.15703 

168.74 

.0637 

13.0 


9.0 


D,/V' 


5.0 


1.0 


0.16 


T/T* 


Figure  5.5.   Dimensionless  hard  sphere  diameter  and 
perturbation  functions  based  on  methane 
reference . 
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nonlinear  regression  routine  and  the  program  GROUPFIT  listed 
in  Appendix  6.   However,  the  data  comparisons  will  normally 
be  made  in  terms  of  the  quantity  AP.   This  procedure  is 
common  and  is  used  to  eliminate  bias  from  the  error  mini- 
mization. 

The  first  test  performed  was  to  calculate  Ap  for  argon 
using  the  characteristic  parameters  as  the  critical 
parameters  and  the  functions  given  by  equations  5-13  and 
5-14.   The  average  error  was  found  to  be  1.42%.   This  is 
larger  than  expected,  but  if  only  the  liquid  phase  data 
are  considered,  the  error  is  just  0.629%,  and  if  the  isotherm 
that  is  within  10%  of  the  critical  temperature  is  neglected, 
the  error  for  the  liquid  isotherms  is  only  0.34%.   This 
is  within  the  accuracy  of  the  equation  of  state  for  this 
region  of  the  phase  diagram.   On  the  whole  the  largest 
errors  are  found  in  the  low  density  data  near  the  critical 
temperature.   The  equation  is  able  to  reproduce  the  highest 
temperature  isotherms  (T/T   =2.8  and  3.0)  with  an  average 
error  of  0.47%.   So  the  fit  of  the  argon  data  seems  satis- 
factory for  the  desired  usage. 

In  an  attempt  to  improve  the  argon  representation, 
a  regression  was  performed  to  find  the  set  of  characteristic 
parameters  that  minimized  the  sum  square  error  in  the  pres- 
sure changes  as  per  equation  5-17.   The  parameters  found 
were  T*  =  144.45  K  and  V*  =  74.31  cm  /gmole  as  compared 

to  T   =  150.86  K  and  V   =  74.57  cm  /gmole.   The  average 
f  c 
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error  for  all  of  the  data  increased  to  1.53%.   This  worsening 

of  the  percentage  error  occurred  because  the  fit  was  made 

5 
more  even;  that  is,  SSEAp  was  decreased  from  4.549  x  10 

to  3.034  x  103,  but  on  the  average  the  error  was  larger. 
This  is  due  to  the  larger  localization  of  error  in  the 
critical  region. 

A  calculation  was  also  performed  on  the  methane  data 
using  the  argon  functions  and  methane's  critical  parameters 
as  the  characteristic  parameters;  the  error  was  3.13%. 
And  when  a  regression  was  performed  on  the  methane  data, 
the  characteristic  parameters  were  found  to  be  T*  = 
190.54  K  and  V*  =  98.40  cm  /gmole  (compare  Tc  =  190.53 

K  and  V   =98.52  cm3/gmole)  with  an  error  of  2.98%.   This 
c 

shows  that  the  trends  found  for  the  methane  DCFI  as  compared 
to  argon  carries  over  to  the  pressure  change  calculations. 

Similar  calculations  to  those  described  above  were 
performed  using  the  hard  sphere  diameter  and  perturbation 
term  functions  as  given  by  equations  5-15.   The  average 
error  in  pressure  changes  for  all  of  the  methane  data  was 
0.57%  with  the  largest  error  from  the  isotherm  at  175  K 
(within  9%  of  the  critical  temperature).   Even  at  this 
temperature  the  high  density  data  are  well  described. 
In  the  reduced  temperature  range  of  0.5  to  0.89  the  average 
error  was  only  0.18%,  well  within  the  accuracy  of  the  data. 
As  for  argon,  when  the  GROUPFIT  program  was  employed  to 
find  optimal  characteristic  parameters  for  methane,  the 
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SSEAP  was  reduced  from  4.13  x  102  to  2.44  x  10  ,  but  the 
average  error  was  increased  to  0.727%.   Again  .the  localiza- 
tion of  the  error  causes  this  phenomenon. 

A  calculation  of  the  pressure  changes  for  the  argon 
data  using  the  methane  functions  gave  an  error  of  4.65%, 
and  when  optimized,  the  argon  characteristic  temperature 
only  changed  1.3%  and  the  error  reduction  was  only  to  4.56%. 
This  again  shows  that  there  are  noticeable  differences 
in  the  compressive  behavior  of  argon  and  methane  at  high 

densities . 

In  summary,  the  temperature  dependence  of  the  hard 
sphere  diameter  and  perturbation  term  in  the  DCFI  model 
are  seen  to  adequately  represent  the  data  from  which  they 
were  developed.   The  difference  in  DCFI  noticed  for  the 
atomic  species  argon  and  the  polyatomic  methane  are  notice- 
able in  calculations  of  compressions. 

The  rest  of  this  work  will  deal  with  calculations 
of  hydrocarbon  compressions  using  the  group  contribution 
formulation.   The  groups  of  interest  are  methyl  (-CH3) 
and  methylene  (-CH2).   These  are  physically  more  similar 
to  methane  than  to  argon.   For  this  reason  all  of  the  group 
contribution  calculations  will  be  performed  using  the  methane 
reference  functions. 


CHAPTER  6 

USE  OF  THE  GROUP  CONTRIBUTION  MODEL  FOR 

PURE  FLUIDS 


The  development  thus  far  has  been  limited  to  the  formal 
expressions  for  calculating  property  changes  of  mixtures 
based  on  the  group  contribution  model  of  direct  correlation 
function  integrals  and  the  parameterization  of  the  model 
using  data  for  argon  and  methane.   In  both  cases  the  mole- 
cules were  assumed  to  consist  of  only  one  type  of  group. 
This  chapter  presents  the  calculations  of  pressure  changes 
for  pure  fluids  composed  of  different  types  of  groups. 

The  first  section  discusses  the  extension  of  the  cor- 
relations required  for  systems  of  several  groups.   The 
remainder  of  the  chapter  presents  and  discusses  the  results 
of  the  calculations  performed  for  several  n-alkanes  and 
one  n-alkanol. 

Extension  of  the  Model  to  Multigroup  Systems 

The  expressions  developed  in  Chapter  4  for  calculation 
of  pressure  changes  involve  two  contributions. 


,  „HS    A„n PERTURBATION  ,,  , > 

APB  =  AP6    +  AP6  (6-1) 
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The  hard  sphere  contribution  is  calculated  for  the  solution 
of  groups  using  the  results  shown  in  Appendix  5.   To  cal- 
culate the  perturbation  contributions  for  multigroup  systems, 
the  form  of  the  ^  g  functions  must  be  established  for  the 
terms  with  a  ^  B.   In  Chapter  5  a  corresponding  states 
expression  was  developed 


(fa    =  V*   f  (T/T*  )  (6-2) 

cm     aa      '  act 


where  f  is  a  universal  function  of  the  reduced  temperature. 
This  result  is  extended  to  the  unlike  terms  as 

^6  =  Va*B  f  (T/TaV  (6_3) 

and  the  characteristic  parameters  are  found  from 


v*      =  1    ( (V*  )1/3  +  (V*  )1/3)  (6-4a) 

aB    8  l   eta  B6     ; 


T*   =  (T*   T*  )     ( 1-k  „)  (6-4b 

aB      aa  1&B'  aB 


where  k    is  an  empirically  determined  binary  interaction 

aB 

parameter  for  groups  a  and  B.  It  is  important  to  note 
that  this  binary  parameter  can  still  be  determined  from 
analysis  of  pure  component  volumetric  data. 
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Designation  of  Groups 

Before  the  group  contribution  model  can  be  used  for 
calculation  of  property  changes,  the  decision  as  to  which 
collections  of  atoms  in  a  molecule  comprise  a  group  must 
be  made.   The  simplest  choice  would  be  to  consider  each 
atom  as  a  group,  and  thus  the  properties  of  all  molecules 
could  be  described  with  minimal  information.   However, 
that  is  an  overly  ambitious  approach.   Even  for  the  molecules 
considered  here  that  contain  only  carbon,  hydrogen,  and 
oxygen  atoms,  this  is  unworkable.   Typically,  the  common 
organic  radicals  are  designated  as  the  groups. 

For  this  work  it  must  be  decided  if  the  n-alkanes 
are  composed  of  only  one  type  of  group,  or  two.   If  they 
are  considered  as  being  comprised  of  only  one  type  of  group, 
then  the  model  could  never  predict  excess  volumes  for 
n-alkane  mixtures.   This  is  reasonable  at  low  pressures, 
but  for  highly  compressed  systems,  noticeable  excess  volumes 
do  exist  (Snyder,  Benson,  Huang,  and  Winnick,  1974).   Thus, 
the  n-alkanes  will  be  considered  to  be  composed  of  two 
groups,  methyl  (CH3)  and  methylene  (CH2).   To  analyze 
methanol  the  hydroxyl  moiety  (OH)  will  also  be  considered 
as  a  separate  group. 

Pressure  Change  Calculations 

All  of  the  molecules  analyzed  in  this  work  will  be 
considered  to  be  composed  of  only  two  types  of  groups. 
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The  equations  required  for  the  calculation  of  the  pressure 
changes  are 

PB  _  (PB)ref  =  p6HS  .  (PBHS)ref 

-|    (p2-(p2)ref  )(v^11  +  2v1v2^12+V2^22)  (6"5) 

where   here,    v.    is    the   number   of    groups   of    type    i    in   the 
molecule.      The    individual    terms    are    found   using 

pBHS    .    1    , Jo,    +    iilX    +  "* 


"r     'L'*o         (1-C3)2         (I-C3)3 

-    4.2   -^-2- T    }  (6-6a 

(I-C3) 


with 


6-6b 


^9    =    6    P  [vl°l  V2°2     J 


and    the   0 .    terms    are    found   using 


1   03/V*    =    f(T/T*)  (6_6c 

611  1 


where  the  function  of  reduced  temperature  required  in  equation 
6-6c  is  given  in  Chapter  5.   The  perturbation  contributions 
are  found  using  equations  6-3  and  6-4  with  the  required 
function  of  reduced  temperature  again  given  in  Chapter  5. 
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These  functional  forms  are  dependent  only  on  the  reference 
component  chosen  for  their  development,  either  argon  or 
methane.   For  the  two  group  molecules  there  are  then  only 
five  independent  parameters,  V*,  V*,  T*,  T* ,  and  k.^- 

Analysis  of  n-Alkane  Compressions 

In  this  section  we  consider  the  ability  of  the  proposed 
formulation  to  calculate  changes  in  pressure  during  compres- 
sion of  four  n-alkane  molecules:   ethane,  propane, 
n-pentadecane,  and  n-octadecane .   Calculation  of  pressure 
changes  using  the  equations  given  above  were  performed 
for  both  the  argon-  and  methane-based  reference  functions 
as  derived  in  Chapter  5. 

The  initial  plan  of  attack  was  to  determine  the  methyl 
parameters  using  ethane  data  and  ethylene  and  cross  parame- 
ters from  the  propane  data.   The  model  capabilities  would 
then  be  examined  by  using  these  same  parameters  for  calcula- 
tions on  the  long  chain  molecules.   The  analysis  was  per- 
formed both  with  and  without  a  binary  interaction  coeffi- 
cient.  The  characteristic  parameters  for  the  methyl  and 
methylene  groups  as  determined  by  regression  analysis  are 
shown  in  Table  6-1  along  with  the  errors  in  calculated 
pressure  changes.   Several  conclusions  can  be  drawn  from 
these  results  that  will  be  of  future  interest.   First, 
the  best  fit  of  all  of  the  data  was  obtained  using  the 
methane-based  correlations  with  the  inclusion  of  a  binary 
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interaction  coefficient.   This  binary  constant  was  much 
larger  than  usual  and  thus  gave  a  large  value  of  the  T*2 
parameter.   This  results  from  the  fact  that  the  model  is 
much  more  sensitive  to  the  values  of  the  characteristic 
volumes  than  the  characteristic  temperatures  so  that  large 
adjustments  to  the  binary  constant  are  required  to  signifi- 
cantly affect  the  calculations. 

It  was  also  seen  that  2V*.^    =  0.9  Vc>ethane  and 

2v*    +  v*    =  0.9  V  .   This  suggests  an  approach 

CH3     CH2         c, propane 

for  estimating  characteristic  volumes  of  the  groups  from 
critical  properties.   The  characteristic  temperatures  might 
also  be  obtainable  using  critical  values  with  some  sort 
of  mixing  rule  applied. 

While  the  results  are  not  as  good  as  the  molecular 
fit,  they  could  be  satisfactory.   A  more  important  question 
is  their  ability  to  predict  compression  data  for  other 
hydrocarbons.   For  n-pentadecane  and  n-octadecane  results 
were  rather  poor,  the  best  agreement  being  an  average  error 
in  the  pressure  changes  of  16.6%  using  parameter  set  1. 
This  is  discouraging  but  not  unexpected.   The  hard  sphere 
contribution  does  not  behave  properly  at  high  densities 
as  explained  in  Appendix  5.   In  addition,  the  improper 
ideal  gas  limit  of  the  equation  of  state  affects  the  results 

In  an  effort  to  improve  the  model  performance,  the 
data  for  the  two  long  chain  hydrocarbons  were  used  to  deter- 
mine the  methylene  and  interaction  parameters. 
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The  results  of  interest  are  shown  in  Table  6-2.   The 
methyl  parameters  were  those  found  from  analysis  of  the 
ethane  data.   The  results  are  rather  encouraging  in  that 
either  the  four-  or  five-parameter  models  could  give  an 
adequate  representation  of  data  for  both  short-  and  long- 
chain  hydrocarbons.   However,  it  should  be  noted  that  none 
of  the  parameter  sets  could  give  a  reasonable  representation 
of  the  propane  pressure  changes. 

These  results  do  yield  two  interesting  conclusions. 
First,  that  when  the  larger  hydrocarbons  are  considered 
for  the  parameter  estimation  that  the  binary  parameter 
becomes  very  important.   Noticeable  improvement  in  the 
pressure  change  calculations  occurs  when  the  binary  parameter 
is  included.   This  suggests  that  the  methyl  and  methylene 
groups  should  be  considered  as  distinct  entities.   However, 
the  magnitude  of  the  binary  parameter  is  much  greater  than 
would  be  expected.   It  can  also  be  seen  that  the  methane- 
based  correlations  are  more  apropos  for  modeling  the  alkane 
compounds  because  the  groups  are  physically  more  akin  to 
methane  than  to  argon. 

After  studying  the  aforementioned  results,  it  was 
decided  to  establish  a  final  set  of  model  parameters  for 
the  methyl  and  methylene  groups  from  an  analysis  of  all 
of  the  n-alkane  data  simultaneously.   The  methane-based 
reference  functions  were  used,  and  the  calculations  were 
performed  both  with  and  without  a  binary  parameter.   These 
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results  are  shown  in  Table  6-3.   It  can  be  seen  that  use 
of  a  binary  parameter  does  improve  the  results  and  that 
the  order  of  magnitude  of  this  parameter  is  reasonable. 
Using  only  five  adjustable  parameters,  the  average  absolute 
error  in  pressure  changes  calculated  for  all  four  hydro- 
carbons was  5.9%. 

Analysis  of  Methanol  Compressions 

A  set  of  compression  measurements  for  methanol  were 

analyzed  to  determine  if  the  present  model  could  adequately 

calculate  pressure  changes  for  molecules  other  than 

n-alkanes.   Methanol  was  considered  to  be  composed  of  one 

methyl  group  and  one  hydroxyl  (OH)  group.   The  methyl  group 

characteristic  parameters  used  are  those  listed  under 

parameter  set  11  in  Table  6-3.   The  regression  results 

are  shown  in  Table  6-4.   The  characteristic  volume  for 

the  hydroxyl  group  appears  to  be  quite  reasonable  because 

v*    +  v*   =  0  9  V      ,    ,  which  is  the  same  pattern 
CH      OH        c, methanol 

seen  for  the  alkanes.   The  characteristic  temperature  and 
binary  parameter  found  are  surprising,  but  these  values 
were  necessary  to  obtain  a  fit  of  the  data  about  as  good 
as  that  for  the  alkanes. 

Summary 

While  the  calculations  presented  here  are  based  on 
recognized  assumptions  that  are  believed  to  be  sound,  there 
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are  certain  disconcerting  aspects.   In  particular,  the 
equation  of  state  developed  does  not  properly  approach 
the  ideal  gas  limit.   Even  though  the  calculations  are 
always  performed  for  pressure  differences  from  a  liquid, 
or  dense  fluid,  reference  state  this  may  be  important. 

Also,  the  percentage  errors  in  the  pressure  change 
calculations  appear  to  be  rather  large.   In  most  process 
design  calculations  the  temperature  and  pressure  are  known 
and  one  wishes  to  calculate  the  volume.   It  would  then 
be  of  interest  to  know  how  the  present  model  would  perform 
in  calculating  volume  changes  for  given  pressure  changes. 
A  simple  error  analysis  shows  how  this  can  be  estimated. 
At  constant  temperature  errors  in  pressure  calculations 
and  volume  calculations  are  related  by 


AV  =  (f|)   AP  (6-7) 


Then  average  absolute  errors  in  volumes  and  pressures  must 
be  related  by 


/£V/  =  PB  M/  (6-8) 


where  6   is  the  isothermal  compressibility.   For  liquids 
the  product  P6   is  almost  always  less  than  0.2  and  even 
then  is  only  that  large  at  very  high  pressures.   This  sug- 
gests that  if  the  present  correlation  were  used  to  calculate 
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a  volume  change  for  a  given  pressure  change  that  the  average 
absolute  errors  would  be  less  than  20%  of  the  errors  reported 
here  for  the  pressure  changes.   This  shows  that  the  present 
model  compares  reasonably  with  existing  correlations. 

One  reassuring  result  of  regression  analysis  is  the 
magnitude  of  the  group  characteristic  volumes  that  were 
calculated.   Ratios  of  the  group  characteristic  volumes 
found  here  are  very  nearly  equal  to  the  ratios  of  the  van 
der  Waals  volumes  for  these  groups.   This  may  allow  for 
estimation  of  these  group  parameters.   This  type  of  regu- 
larity was  not  seen  for  the  group  characteristic  tempera- 
tures.  Also,  whereas  the  characteristic  volumes  for  the 
different  groups  did  not  vary  much  when  determined  from 
different  data  sets,  the  characteristic   temperatures  did. 
As  a  result,  the  binary  constants  did  not  follow  a  pattern. 
The  binary  constant  was  necessary  in  all  cases  to  obtain 
an  optimal  representation  of  the  data,  but  the  improvement 
obtained  by  introducing  this  additional  parameter  was 

minimal. 

No  calculations  have  been  reported  for  chemical  poten- 
tial changes  or  for  pressure  changes  for  multicomponent 
systems.   The  necessary  formulae  have  been  presented  in 
Chapter  4  and  Appendix  5.   Pressure  change  calculations 
for  multicomponent  systems  are  no  more  difficult  than  those 
reported  here  because  the  molecules  considered  contained 
more  than  one  group. 


CHAPTER  7 
DISCUSSION 


This  work  has  centered  on  the  development  of  a  group 
contribution  liquid  phase  equation  of  state.   Fluctuation 
solution  theory  was  used  to  develop  the  exact  relationships 
between  thermodynamic  derivatives  and  direct  correlation 
function  integrals.   Two  approximations  were  made  to  actually 
construct  the  equation  of  state.   The  first  approximation 
made  was  the  use  of  the  interaction  site  formalism,  or 
RISM  theory.   The  last  stage  in  the  model  formulation  was 
the  choice  of  an  approximate  form  for  the  group  direct 
correlation  function  integrals.   After  that  development, 
the  general  expressions  for  calculation  of  pressure  and 
chemical  potential  changes  were  written  and  used  for 
several  compounds. 

In  this  chapter,  the  important  aspects  of  each  stage 
in  the  development  described  above  will  be  examined  in 
detail.   Great  emphasis  will  be  placed  on  the  approximations 
embodied  in  the  model  development  as  they  pertain  to  the 
numerical  results. 
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Fluctuation  Solution  Theory 

The  relationship  between  integrals  of  molecular  correla- 
tion functions  and  thermodynamic  derivatives  have  come 
to  be  known  as  fluctuation  solution  theory  (O'Connell, 
1981).   Kirkwood  and  Buff  (1951)  originally  reported  these 
relationships  though  they  have  not  been  extensively  used. 
The  original  results,  using  total  correlation  function 
integrals,  are  valid  for  both  spherical  and  molecular  sys- 
tems.  With  the  introduction  of  the  direct  correlation 
functions,  through  the  Ornstein-Zernike  equation,  the  analy- 
sis is  not  as  simple.   Appendix  1  details  a  procedure  which 
surmounts  all  difficulties  encountered.   This  result  is 
not  new,  but  the  rigorous  proof  appears  to  be.   This  result 
could  also  be  derived  through  the  use  of  the  operator  tech- 
nique of  Adelman  and  Deutch  (1975). 

Approximations 

While  equations  2-24  and  2-25  are  exact  results  for 
fluctuation  derivatives,  one  must  have  values  for  the  direct 
correlation  function  integrals  for  these  to  be  useful. 
Because  exact  values  of  the  required  integrals  are  not 
generally  available,  some  approximations  must  be  used  to 
obtain  these  values.   The  approximations  used  will  determine, 
to  a  large  degree,  the  success  of  any  thermodynamic  calcula- 
tions.  In  this  work,  two  levels  of  approximation  are  used 
and  they  will  be  discussed  separately  below. 
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Group  Contributions 

One  of  the  major  goals  of  this  research  was  the  develop- 
ment of  a  group  contribution  formulation  of  fluctuation 
solution  theory.   This  appears  to  be  a  realistic  and  rational 
goal.   This  is  true  because  fluctuation  solution  theory 
requires  knowledge  of  correlation  functions  and  correlations 
among  the  groups  are  well  defined.   The  major  problem 
encountered  was  the  ability  to  separate  the  inter-  and 
intramolecular  correlations.   The  RISM  theory  (Chandler 
and  Andersen,  1972)  purports  to  affect  this  separation 
and  offer  an  approximate  relationship  between  group  direct 
correlation  functions  and  molecular  correlations.   The 
calculations  reported  here  are  based  on  this  original  version 
of  the  RISM  theory. 

Recently,  the  direct  correlation  functions  calculated 
using  the  RISM  theory  have  become  the  objects  of  considerable 
attention  (Cummings  and  Stell,  1981,  1982b;  Sullivan  and 
Gray,  1981).   It  has  been  shown  that  in  many  situations 
for  neutral  molecules  that  the  site-site  direct  correlation 
functions  are  necessarily  long-ranged.   This  means  that 
these  functions  are  nonintegrable  and  thus  no  compressibility 
relation  involving  these  functions  may  exist.   Cummings 
and  Stell  (1982a)  have  shown  how  a  general  compressibility 
relation  can  be  derived  in  terms  of  a  limiting  operation, 
but  this  is  not  a  practical  solution.   However,  for  the 
cases  that  have  been  considered  in  the  literature  it  seem 
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this  problem  can  be  overcome.   If  the  compressibility  rela- 
tion is  viewed  as  the  k  ->  o  limit  of  a  result  in  Fourier 
space,  where  all  of  the  correlation  function  transforms 
exist,  then  all  of  the  mathematical  manipulations  required 
may  be  performed.   Then  as  the  long  wavelength  limit  is 
taken,  we  find  (Appendix  7)  the  following  relationship: 

f3pe1      =  L      yyx.x.yy  v.  «.fl  c  fl(o  (7-d 

l3p  J^      p  U    i  D^   la  DB   aB 

While  the  individual  functions  Cag(r)  may  be  long-ranged, 
it  appears  that  the  collection  of  terms  u>cw  is  well  defined 
in  the  k  •*■   o  limit.   It  is  known  (Cummings  and  Stell,  1982a) 
that  for  diatomic  molecules  this  is  true,  and  thus  a  simple 
compressibility  theorem  applies.   Also  for  triatomics, 
as  examined  by  Cummings  and  Stell  (1981),  we  find  that 
this  summation  of  terms  is  nondivergent .   As  has  now  become 
known,  the  possible  divergencies  in  the  Cag(r)  functions 
are  due  to  intramolecular  effects  that  are  not  shown 
explicitly.  It  appears  that  the  projection  by  the  v  matrix 
removes  the  divergent  terms.   This  is  analogous  to  the 
ionic  solution  case  (Perry  and  O'Connell,  1984)  where  the 
charge  neutrality  constraint  removes  the  divergent  part 
of  the  direct  correlations. 

Chandler  et  al.  (1982)  formulated  the  proper  integral 
equation  theory  for  site-site  correlations  to  offer  an 
exact  formulation  involving  nondivergent  Cag(r)  functions 
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that  here  are  labeled  C°g(r).   These  new  functions  are 
related  to  the  RISM  C   (r)  functions  and  are  made  nondiver- 
gent  by  removal  of  intramolecular  correlations  involving 
only  sites  a  or  6.   This  requires  use  of  auxiliary  functions 

ft  „(r)  that  have  a  complicated  density  dependence.   These 

aB 

functions  do  possess  the  interesting  property 

ftaS(o)  =  ft  ¥a,6  (7-2) 

It  is  then  possible  to  derive  a  compressibility  relation 

involvinq  only  C°  (r)  function  integrals  and  ft 

aB 


3PS 


i      =  k  -  p  yyx.x.yyv •  v .D 

T,X  i]    JaB     J 


J  d?  C°   (r)  (7-3 


In  this  case,  the  integrals  always  exist.  However,  we 
have  no  knowledge  of  the  functional  dependence  of  ft  on 
temperature  or  density. 

One  further  aspect  of  the  group  contribution  approach 
used  here  needs  to  be  mentioned.   This  formulation  works 
on  the  Ansatz  that 


c^d.2)  =  11    ^iavj6ca6(?a6) 

J        a  B 


7-4 
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which  is  only  to  be  considered  an  approximate  relation. 
However,  the  compressibility  relation  obtained  from  the 
above  starting  point  and  that  found  using  the  proper  integral 
equation  formalism  (as  shown  in  Appendix  7)  are  the  same. 
The  only  difference  is  that  all  forms  of  the  RISM  theory 
assume  that  each  group,  as  opposed  to  each  type  of  group, 
have  separate  correlation,  i.e.,  that  the  set  of  functions 
ca?(r)  should  be  considered,  not  just  cag(r).   This  is 
definitely  true  at  the  level  of  the  correlation  functions 
themselves.   However,  it  is  always  possible  to  define  a 
set  of  c  0(r)  functions  that  are  averages  of  the  c. .(r) 

CLP  XJ 

functions  (Adelman  and  Deutch,  1975)  and  retain  the  same 
compressibility  relations. 

Modeling  of  Direct  Correlation  Function  Integrals  (DCFI) 

As  the  DCFI  are  functions  of  temperature  and  density, 
it  would  be  possible  to  model  them  using  polynomial  expan- 
sions for  the  pure  fluids  and  a  solution  theory  for  the 
mixture  quantities.   In  this  work,  we  attempted  to  use 
f luctuational  forms  that  have  a  theoretical  basis.   This 
approach  seems  to  have  both  merits  and  demerits.   For  mole- 
cules, Mathias  (1979)  had  success  in  modeling  DCFI  using 
an  approach  suggested  by  perturbation  theory.   In  the  present 
case  this  is  not  of  the  greatest  value  because  analytic 
forms  are  available  for  group  DCFI  in  only  some  limited 
cases  (Morris  and  Per  ram,  1980 ;  Cummings  and  Stell,  1982a), 
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and  even  then,  the  solutions  are  not  in  closed  form.   Thus, 
some  reasonable  assumptions  must  be  made  to  proceed  with 
the  modeling. 

It  is  known  that  the  RISM  correlation  functions  can 
be  related  to  Percus-Yevick  correlation  functions  for 
molecular  species  in  some  limiting  cases  (Chandler,  1976). 
This  led  to  the  use  of  the  form  for  the  DCFI  presented 
in  Chapter  5.   The  important  aspect  of  this  approach  is 
that  the  groups  are  dealt  with  as  if  they  were  all  indepen- 
dent.  This  is  not  physically  the  situation.   The  dependence 
is  caused  by  the  intramolecular  correlations,  and  they 
must  be  properly  accounted  for  to  yield  an  accurate  model. 
Originally  it  was  believed  that  the  RISM  theory  took  account 
of  this  behavior,  but  as  was  mentioned  in  the  preceding 
section,  this  is  not  true.   Even  the  proper  integral  equation 
direct  correlation  functions,  equation  7-3,  while  short 
ranged,  still  involve  intramolecular  correlations  due  to 
third  body  effects.   Cummings  and  Sullivan  (1982a)  explicitly 
show  that  while  the  hag(r)  functions  are  purely  intermolecu- 

o   "*" 

lar  quantities  that  the  cag(r)  are  not. 

The  exact  resummation  of  the  cluster  expansion  for 
h  g(r)  shown  by  Chandler  et  al.  (1982)  suggests  a  way  around 
the  aforementioned  difficulties.   It  should  be  possible 
to  define  a  set  of  direct  correlation  functions,  labeled 
here  as  c  g(r),  that  include  only  intermolecular  effects. 
Chandler's  (1976)  analysis  can  be  examined  to  see  that 
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these  diagrams  can  be  isolated.   And  then,  in  the  same 
fashion  as  was  used  to  eliminate  the  long-ranged  behavior 
of  the  original  RISM  direct  correlation  functions,  it  would 
be  possible  to  write  an  exact  proper  integral  equation 
of  the  form 


w  +  Ph  =  [I  -  PWc  ]   W  ( 7-5> 


where  the  elements  of  the  W  matrix  would  contain  all  infor- 
mation about  intramolecular  effects.   It  is  then  easy  to 
show  the  relationship  between  the  different  classes  of 
direct  correlation  functions,  for  example 

g  =  i  +  -  (i_1  -  i"1)  (7_6) 

The  elements  of  the  W  matrix  would  be  functions  of 
temperature  and  density,  but  the  exact  functionality  could 
not  be  determined  in  general.   It  is  possible  to  identify 
some  of  the  properties  of  the  new  functions,  such  as 


piomW^  (7-7) 


and 


,  lim  *      T^r   tt  <  7-8 ) 

k->-n   aJo=W~Vn  k  /     o  > 

a6     aB   a ,  8 
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Using  equations  7-5  and  7-8,  in  conjunction  with  the 
techniques  of  Appendix  7,  a  compressibility  relation  can 
be  derived  in  terms  of  these  new  group  direct  correlation 
functions 


3PB 
3p 


i       =  -  -  p  yyx.x.yyv •  v  -a 

T,X  l]    JaB     J 


J  d?  C \    (r)  (7-9 


All  of  the  above  analyses  show  that  there  is  no  simple 
solution  to  the  modeling  problem.   One  may  either  work 
with  a  simple  compressibility  relation  written  in  terms 
of  DCFI  whose  behavior  is  not  certain,  or  in  terms  of  well 
behaved  DCFI  but  have  another  unknown  quantity  present. 
The  former  approach  was  taken  for  this  work  and  is  discussed 
below. 

DCFI  Model 

The  rationale  behind  the  development  of  the  DCFI  model 
has  been  presented  in  Chapter  5.   It  is  of  interest  to 
note  that  the  proposed  form  is  a  group  variation  of  a  van 
der  Waals,  or  mean  field,  model.   Appendix  2  details  how 
higher  order  correction  terms  could  be  appended  to  the 
proposed  form. 

Because  the  DCFI  model  was  formulated  as  a  correspond- 
ing states  correlation,  the  constants  in  the  temperature 
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dependence  of  the  hard  sphere  diameters  and  perturbation 
terms  had  to  be  determined  using  experimental  data.   Volu- 
metric data  for  both  argon  and  methane  were  used  to  accom- 
plish this  task.   The  results  of  the  compression  calculations 
show  that  the  methane-based  correlations  were  definitely 
superior  for  representation  of  the  properties  of  the  larger 
molecules . 

Comparison  Calculation 

Chapter  6  presents  the  results  of  the  group  contribution 
compression  calculations  for  the  n-alkane  and  methanol. 
The  results  seem  to  indicate  that  the  theoretical  basis 
is  feasible  but  that  the  accuracy  is  not  that  desired. 
It  is  an  accomplishment  to  be  able  to  perform  the  compression 
calculations  for  all  of  the  n-alkanes  using  only  five 
parameters,  but  the  overall  accuracy  is  not  high.   One 
would  hope  to  be  able  to  calculate  pressure  changes  with 
an  accuracy  of  about  1%  (this  allows  for  density  change 
calculations  accurate  to  -0.2%),  and  this  cannot  be  assured. 

The  theoretical  basis  behind  the  equation  of  state 
development  is  now  on  firm  footing.   Our  compressibility 
relation  is  exact;  one  must  simply  have  proper  models  for 
the  functions  involved.   Here,  this  appears  to  not  be  the 
case.   Because  we  use  the  RISM  approximation,  the  cQg(r) 
functions  must  contain  intramolecular  effects  that  have 
no  analogy  at  the  molecular  level.   Thus,  our  model  that 
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is  based  on  an  analogy  to  a  previously  successful  molecular 
approach  appears  inadequate.   Note,  however,  that  we  seem 
to  have  little  recourse.   Because  if  we  were  to  introduce 
easily  modelable  c0g(r)'s,  then  the  aspect  of  determining 
the  W  function  must  be  addressed,  and  this  is  still  an 
unsolved  problem. 

Summary 

This  study  has  revealed  many  interesting  aspects  about 
site-site  correlation  functions  and  their  relationship 
to  thermodynamic  derivatives.   A  group  contribution  liquid 
reference  state  equation  of  state  has  been  formulated  and 
tested.   The  results  are  encouraging  enough  to  warrant 
further  investigation  but  not  accurate  enough  for  practical 
density  calculations.   It  appears  that  the  major  problem 
associated  with  the  present  approach  is  the  inability  to 
model  the  required  functions,  and  suggestions  have  been 
made  as  to  how  one  can  alleviate  that  problem. 


CHAPTER  8 
CONCLUSIONS 


The  goal  of  this  work  was  the  development  of  a  group 
contribution  technique  for  calculation  of  liquid  phase 
properties  of  mixtures.   The  general  results  are  available 
but  have  only  been  tested  on  a  few  pure  components.   The 
calculated  pressure  changes  for  four  n-alkanes  and  methanol 
are  in  general  in  error  by  less  than  7%.   This  is  reasonable 
but  not  of  high  enough  accuracy  for  process  engineering 
calculations.   While  the  results  of  this  study  are  not 
outstanding,  there  are  several  interesting  conclusions 
that  can  be  drawn. 

It  has  been  shown  that  it  is  possible  to  derive  thermo- 
dynamic property  derivatives  from  a  RISM  theory.   A  general- 
ized compressibility  theorem  was  proven,  and  this  was  used 
along  with  a  model,  based  on  a  rigorous  perturbation  theory, 
to  develop  a  van  der  Waals  equation  of  state. 

The  liquid  phase  equation  of  state  was  applied  to 
the  representation  of  compression  measurements  of  argon 
and  methane.   It  was  shown  that  proper  choice  of  the 
repulsive  contribution  to  the    equation  of  state  is 
important.   Also,  that  an  optimum  repulsive  contribution 
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could  be  determined  from  a  new  hard  sphere  equation  of 

state. 

The  use  of  the  group  contribution  model  was  limited. 
Pressure  changes  could  be  calculated  for  both  long-  and 
short-chain  alkanes  using  a  five-parameter  model.   The 
results  were  of  reasonable  accuracy  for  these  cases  and 
when  the  model  was  applied  to  methanol.   While  the  volumetric 
parameters  in  the  model  appeared  to  correlate  well  with 
previous  molecular  results  (Mathias,  1979),  the  temperature 
parameters  followed  no  pattern. 

An  analysis  of  the  present  work  also  shows  several 
areas  that  require  future  work.   First,  the  present  model 
should  be  applied  for  calculation  of  chemical  potential 
changes  for  n-alkane  mixtures.   This,  along  with  mixture 
volumetric   calculations,  would  be  a  strong  test  of  the 
model's  ability.   Secondly,  some  effort  to  identify  the 
intramolecular  correlation  functions  of  the  form  shown 
in  Chapter  7  may  be  required  to  improve  the  present  model. 
This  appears  to  be  the  weakest  aspect  of  the  present  work. 


APPENDIX  1 

FLUCTUATION  DERIVATIVES  IN  TERMS  OF  DIRECT  CORRELATION 

FUNCTION  INTEGRALS 

To  eliminate  the  difficulties  associated  with  the 
form  of  the  general  0-Z  equation  we  define  an  angle  averaged 
total  correlation  function,  <hi.(R-L,R2)>  by 

<h.  .(R,,R„)>  =  -\    j    dfi,dfi?  h   (1,2)  (Al-1) 

13   1   2      ^2      1   2   13 

and  then  we  implicitly  define  a  new  set  of  direct  correlation 
functions,  the  <c . .(R1,R2)>  by  an  0-Z  type  equation 


<h..(R1,R2)>  =  <ci;j(RrR2)> 


<N,  > 
+  I  -^-  J  dR3<cik(R1,R3)><hkj(R3,R2)>   (Al-2) 

k 


Now,  because  the  full  h.  .  (1,2)  is  translationally  invariant, 
the  average  function  must  also  be,  and  thus  the  <c^±> 
function  must  also  possess  this  property.   Then  if  we 
integrate  equation  Al-2  over  one  coordinate  we  have 


<h.  . (R)>  =  <c.  .  (R)> 

13      <Nv>1 

+  I    — =£-  J  ds  <c,  (s)Xh,  .  (R-s)>         (Al-3 
!r        V   '  lk       k] 

And  if  this  is  then  integrated  once  more  we  find 
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<h.  •>  =  <c  .>  +  Y  <N,  ><c..Xh,  .>  (Al-4) 

lj         lj      L  k     IK     K] 

k 


where  we  define 


<h.  .>  =  |  J  dR  <hi.(R)>  (Al-5a) 


<c.  .>  =  h    J  dR  <c.  .(R)>  (Al-5b! 


and  when  these  terms  are  collected  in  matrix  form  we  find 
<h>  =  <g>  +  <g><g><h>  (Al-6) 

Here,  the  terminology  is  the  same  as  that  used  in  Chapter 
2  for  the  definition  of  the  matrix  elements.   This  relation 
becomes  most  useful  when  it  is  realized  that 

<h>  =  H  <Al-7) 

as  can  be  shown  by  direct  substitution.   And  now  rearrange- 
ment leads  to 


[N  +  NHN]  X  =  [N  1    -    <c>]  (Al-8) 


which  shows  that  the  fluctuation  derivatives  can  be  written 
in  terms  of  the  integrals  of  the  <^i.>  functions.   The 


relation  between  c±  •  and  <cj_.>  is  not  obvious  but  can  be 


■  D        ]  3 
derived.   By  definition 


<c.  .  (R, ,R,)>  =  <h.  . (R, ,R,)> 

ij         1   2        i]   X   z 


-I-Tldl3  <cik(«1,S3><hkj(«3,«2)> 

k 


and,  if  the  0-Z  equation  is  angle  averaged,  we  find 

<hij(R1,R2)>  =  ^  J"  d¥"2  cij(1'2) 

+  I   __*_   J  d3dX?1d^2  cik(l,3)hkj(3,2)         (Al-9) 
k  Vft 

Now,  when  these  two  relations  are  combined,  the  link  between 
the  angle  averaged  c± .  and  the  <ci->  functions  are  given 


<ci.(R1,R2)>  =  \    J  d^d^2  Ci:j(l,2) 

+    Y   —f-    J    d?d^1dJT2    cik(l,3)hk.(3,2) 
k   Vfi 

-    I    ^    MR3    <cik(R1,R3)><hkj($3,R2)>  (Al-10) 

k 

And  now  insert  the  definition  of  <hk^(^3,R2)>  to  find 
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<c,4(R1fR0)>  =  -\    J  d^d^0  c.  .,(1,2) 


<1 

+  I   -f-   J  d3d^,d^9  c.k(l,3)h   (3,2) 
v  „nJ  1   ^   IK        KJ 


k  vsr 


-  I  ^T"  /  d«2<cik(Rl,R3)>hkj(3,2) 


and  the  terms  can  be  grouped  as 


(Al-11) 


<c.  .(R,,lL)>  -    J  d^,d^„  c.  .(1,2) 
13   1   2      n      1   2   i] 

<N,  >     _^  ,  ^   + 

=  I    -~-    J  d3dJT2hk,(3,2)  {£  Jd^1  cik(l,3)-<cik(R1,R3 
k  Vfi  J 


)>} 


(Al-12) 


Now,  use  the  fact  that 


:ci.(RrR2)>  ^!dfi2  <c   (RrR2)>  (Al-13 


to  rewrite  equation  Al-12  as 

I  /  d^  {<c.j(R1,R2)  -^  J  dS]_  c..(l,2)} 

=  I  j¥  J  «d^  hkj(3,2)  {§  J  d^  cik(l,3 


-  <c.k(R1,R3)>: 


:ai-i4 


It  is  now  advantageous  to  define 


For  now  equation  Al-14  takes  the  form 


S  >  dfi2  Kij(5l'52^2) 
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;.  .(R^,^)  =  <c.  j(R1,R2)>  -  |  J  d\  c.  .(1,2)     (Al-15 


<N,  > 


=  £  — ^S_   /  33d^2  hk .(3,2)Kik(R1R3^3)         Al-16) 


k  vn 


and  through  the  use  of  the  identity 


11   1  z   ^ 


=  [  5jk  J  d3  Kik(R1R3n3)  6(R2-R3)6(^2-i!3) 


(Al-17 


It  is  now  possible  to  write  the  expression  for  K^ .  as 


I    -^  J  d3d^2  K\   (R-jR,^)  6(r\,-R\)  6  (&,-$,  ) 


'3  3'   v  2   3'     2   3 


<Nk> 


=  -  I   — £-  f  d3dfi„  h,  .(3,2)  K.,  (R,R,fl, 

K  VQZ  l       KD 


'Al-ll 


and  this  can  be  rearranged  to 
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I    1  J  d3dfi2  Kik(RrR3^3)   {6jk  6(R2"R3)  6(fi2"fi3) 

k 


<N  > 

+  "W  hkj 


k   -   (3'2)  }  =  0      (Al-19) 


The  term  in  brackets  is  multicomponent  density  fluctuation 
functions,  which  is  nonzero  almost  everywhere  (Evans,  1979 
And  thus  a  sufficient  condition  for  equation  Al-19  to  be 
valid  is  that 


(Al-20) 


(Al-22 


Kik(Rl'R3'V  =  ° 


If  this  is  true,  then  it  is  easy  to  show  that 
<g>  =  g 

which  was  the  desired  result. 

And  thus  equation  Al-8  can  be  written  as 


[N  +  NHN]  l    =    [N  l    -    g]  (Al-23 


APPENDIX  2 
PERTURBATION  THEORY  FOR  DIRECT  CORRELATION  FUNCTION 
INTEGRALS  USING  THE  RISM  THEORY 

In  Chapter  4  it  was  noted  that  theoretical  results 
from  perturbation  theory  can  be  useful  in  determining  func- 
tional forms  for  empirical  models  of  direct  correlation 
function  integrals.   An  exact  result  for  the  direct  correla- 
tion function  perturbations  from  perturbations  in  total 
correlation  functions  based  on  the  RISM  theory  is  derived 
below. 

The  problem  considered  here  is  to  find  the  perturbation 
in  Cag,  designated  as  C^g,  that  is  induced  by  a  change 
in  h^g  from  h°6  to  ha°   +  ha].      We  begin  with  the  Fourier 
transform  of  the  RISM  form  of  the  Ornstein-Zernike  equation 

K,    =  H  Z         c    «    +  PQ  II    %y  cyn  hnB  (A2-1) 

The  caret  O  denotes  the  Fourier  transform  and  the  functions 
u    are  the  transforms  of  intramolecular  correlation  func- 
tions.   We  rewrite  equation  A2-1  in  matrix  form  as 


H  =  ucu  +  p  ^.QU 


Equation  A2-2  can  be  rearranged  to 
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(A2-2) 
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[I  -  p  ug][I  +  P0M    l]    =  I  (A2-3) 

where  I  is  an  identity  matrix.   The  terms  in  the  matrix 
H  are  split  into  reference  contributions  (g°)  and  perturba- 
tion contributions  ( g1 )  where  the  g°  are  assumed  to  obey 
an  equation  like  A2-3  with  a  matrix  of  reference  contribu- 
tions to  the  direct  correlation  functions  (g  ).   We  wish 
to  determine  the  perturbation  induced  in  the  direct  correla- 
tion functions  (c  )  defined  by 


gl  =  £  _  £o  (A2-4) 


Use  of  a  mathematical  identity  leads  to 

[I  -  p  uc]  =  [|  +  P  I°^~1]  ~  [i  +  Po=°=~    ]    Po=  = 


*  [I  "  P  wc]       (A2-5 


Solution    for    c      gives 


IQ1    =    [|    +    (J"p0^£0)po=    =       ^ 


[(I-P    wc°)p    H1^    1(I    -    P^c°)]  (A2-6) 

N=      o==         o=    =         =  o 


This  result  is  formidable,  but  exact.   It  can  be  expanded 
in  powers  of  g   to  yield 
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we,  =  H1^  1  -  p   §  w   uco  -  PQ  ^cj  a 


+  p  2  wc  H1^"1^^   +  0[(H   )  ]  (A2-7) 

o   ==o=    -   o 


This  result  is  still  difficult  to  use.   However,  the  low 
density  limit  is  given  by 


iili  =  pj~°  i1  (A2_8) 


This  is  a  RISM  form  of  the  mean  spherical  approximation. 
In  the  low  density  limit  we  have 


h  fl  =  J  dld2d?6(r?)6(H  -r)  f(l,2) 


e  k'r       (A2-9 


where  f(l,2)  =  C"Bu(1'2)  -  1  and  u(l,2)  is  the  intermolecular 
pair  potential.   If  u  is  split  into  reference  and  perturba- 
tion terms,  we  find  that 


1 


-ik  *r 


J  dld2d?5(?«)  6 <?£-?)  C"6U  (1'2)  f1(l,2)e-lk* 


(A2-10 


If  the  reference  system  is  a  hard  sphere,  or  hard  core, 
system,  then  the  result  in  equation  A2-10  simplifies  to 

h  *  =  /  dld2  6(?°)  6(?J-?J  f1(l,2)e-k'?  (A2-11) 
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If  the  intermolecular  perturbation  potential  is  written 
as  a  sum  of  group-group  terms 

u^l^)  =  H   u\x    ir[-rX2)  (A2-12) 

Then  the  elements  of  the  H  martrix  can  be  found  as 


H1  =  I       (^(-Bu1) )  fi/nl  (A2-13 

n=l   ~ 


This  allows  for  identification  of  the  perturbation  in  the 
direct  correlation  function  as 


c 


_1  _  _  Bu1  +  B2wu1^/2!  -  B^uiu1  w/3  !  +  .  .  .       (A2-14) 

And  then,  at  a  macroscopic  level  (k+o)  this  suggests  that 

C1   =  a    +  S   /T  +  c   /T2  +  .  .  .  (A2-15) 

aB     ay     ay       aB 

The  perturbation  should  be  an  expansion  in  powers  of  1/T. 
This  form  is  utilized  in  this  work. 


APPENDIX  3 
THREE  BODY  DIRECTION  CORRELATION  FUNCTION  INTEGRALS 

It  is  possible  to  define  a  hierarchy  of  direct  correla- 
tion functions  and  their  integrals  for  any  number  of  mole- 
cules.  In  specific,  the  three  body  direct  correlation 
function  integrations  are 


9C. 


C 


ijk    l3pk  ;T,p 


JJh  (A3-1) 


k   x ' H £^k 


If  the  two  body  direct  correlation  function  integral  is 
represented  using  the  RISM  theory 


Cij  =  H  viaVjBCa| 


(A3-1.5) 


Then  the  total  differential  of  this  function,  at  constant 
temperature,  is 

dC.  =  HI    v.  v    (^2!)        d  (A3-2) 

When  a  three  site  direct  correlation  function  integral 
is  defined  in  a  manner  analogous  to  equation  A3-1,  then 
a  relationship  between  the  molecular  and  site  three  body 
functions  is  obtained 
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9p 

c. .,  =  in  v.  v <0v.  c  _  (^)  (A3-3 

ink    LLDL       ia  :8  ky  aB  v3pv'p0 
apy  K   ^ 


and  the  use  of  the  material  balance  then  leads  to 


C  .,  =  [H  v.  v  ..v,  C  a  (A3-4) 

ink    L%L      ia  38  ky  aBy 


This  simple  relationship  shows  how  the  RISM  approxima- 
tion can  be  extended  to  multibody  correlation  functions. 
It  is  obvious  that  this  technique  can  be  extended  ad  nauseam 
to  form  group  contribution  theories  for  higher-order 
correlation  functions. 


APPENDIX  4 
NONSPHERICITY  EFFECTS 


Experimental  data  presented  in  Chapter  5  (that  was 
used  to  develop  figures  5-1  and  5-2)  show  that,  for  dense 
liquids,  a  two-parameter  corresponding  states  relationship 
exists  as 


|PB  =  f(p,  t)  (A4-1 

3p 


This  would  be  expected  for  simple  molecules  with  spherically 
symmetric  potential  functions.   The  fact  that  this  same 
correlation  can  apply  to  simple  molecules  (Ar),  long  chain 
molecules  (normal  hexadecane )  and  polar  species  (ammonia) 
is  somewhat  surprising.   This  appendix  offers  justification 
for  this  type  of  correlation  using  both  microscopic  and 
macroscopic  agreements. 

General  Considerations 

The  intermolecular  potential  determines  the  state 
of  a  system  for  a  given  density  and  temperature  in  the 
absence  of  external  fields.   The  actual  force  that  acts 
on  any  molecule  is  dependent  on  the  phase  space  coordinates 
of  all  other  molecules  in  the  system. 
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For  simplicity  these  initial  considerations  will  be 
discussed  in  terms  of  a  symmetric  pair  potential.   For 
many  molecules  this  is  similar  to  that  shown  in  Figure 
A4-1.   For  intermolecular  separations  less  than  r^    the 
interaction  is  repulsive,  and  the  interaction  is  attractive 
for  greater  separations.   For  dense  fluids  the  typical 
intermolecular  separation  is  less  than  r^    so  that  the  domi- 
nant forces  are  repulsive.   For  low  pressure  gases  the 
typical  separation  is  of  the  magnitude  of  r^  or  larger. 
Then  the  dominant  type  of  interaction  is  attractive. 

This  is  surely  a  simplified  picture  of  reality  because 
for  real  molecules  the  orientations  must  be  considered. 
There  are  both  short-ranged  and  long-ranged  forces  of  the 
attractive  and  repulsive  forms.   However,  this  suggests 
the  proper  form  to  follow  for  the  analysis.   There  are 
two  separate  contributions  to  consider: 

1.  Molecular  shape--the  short  range  repulsive  forces 
should  depend  on  the  shape  of  the  molecule. 

2.  Long-range  anisotropic  f orces--multiple  moments 
are  of  long  range  and  depend  on  molecular  orientation. 

Molecular  Shape  Effects 

For  dense  fluids  the  molecules  are  tightly  packed. 
The  packing  is  determined  by  the  harshly  repulsive  part 
of  intermolecular  potential.   For  nonspherical  molecules 
these  effects  must  be  shape  dependent.   Barker  and  Henderson 
(1981)  and  Steele  and  Sandler  (1974)  have  shown  that  harshly 
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Figure  A4-1.   Typical  intermolecular  potential 
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repulsive  potentials  can  be  modeled  as  hard  core  potentials 
if  the  core  size  is  taken  as  temperature  dependent.   The 
analysis  here  will  then  be  to  compare  hard  sphere  systems 
with  general  hard  body  systems. 

A  general  hard  body  theory  does  not  exist.  The  proper- 
ties of  hard  convex  bodies  have  been  determined  using  the 
scaled  particle  theory  (Gibbons,  1970)   and  an  extension 
of  the  Carnahan-Starling  theory  (Boublik,  1974).   Rigby 
(1976)  has  shown  that  it  is  possible  to  calculate  the  prop- 
erties of  concave  bodies  using  the  convex  body  equations 
if  the  parameters  are  properly  selected. 

Two  points  are  to  be  addressed  here.   First,  can  general 
convex  bodies  be  represented  as  spheres?   An  affirmative 
answer  to  this  question  could  explain,  to  a  large  degree, 
the  existence  of  the  corresponding  states  principle  for 
the  dense  fluids.   A  second  question  to  consider  is  will 
the  proposed  theory,  that  for  overlapping  spheres,  give 
a  better  representation  of  the  molecular  properties  than 
the  spherical  model. 

For  hard  convex  bodies  Boublik  (1974)  has  presented 
what  seems  to  be  the  most  accurate  expression  (Boublik, 
1975).   For  a  pure  component  this  is 

7  _  l+(3Y-2)y+(3Y2-3  +l)y  -y  y  (A4-2) 

(1-y) 
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where 


y  =  PV 

V  =  molecule  volume 

Y  is  a  single  parameter  that  represents  the  non-sphericity. 

For  had  spheres  Y  =  1  and  equation  A4-2  reduces  to  that 

of  Carnahan  and  Starling.   To  determine  the  ability  of 

the  hard  sphere  equation  to  represent  the  hard  convex  bodies, 

for  calculation  of  pressure  changes,  the  problem  is 

formulated  by  realizing  that  in  general  from  equation  A4-2 

we  have 


3pB 
3p 


f (y,Y) 


(A4-3) 


An  equivalent  hard  sphere  size  was  found  that  minimized 
sum  squared  errors  in  pressure  deviations  calculated  over 
the  dense  fluid  range.   The  following  variational  problem 
was  solved: 

Find  a  such  that 


[f (y,Y)  -  f (ay,l)  ]  dp 


A4-4) 


o. 


was  a  minimum.   The  value  of  a  was  found  to  be  sensitive 
to  the  density  range  evaluated  but  the  following  values 
were  used  as  representative. 


P  .   =  0.2/V 

mm       — 


A4-5a 
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P     =  0.6/V  (A4-5b) 

max       —  .  . 

The  value  of  a  was  found  to  be  a  smooth  function  of  y   for 
values  of  Y  from  1  to  3 .   Y  =  3  would  correspond  to  a  sphero- 
cylinder  with  a  length  to  diameter  ratio  of  about  10. 
Figure  A4-2  shows  the  comparison  for  a  spherocylinder  system 
using  these  results.   Notice  that  the  agreement  is  not 
exact.   The  equivalent  hard  sphere  is  always  one  with  a 
volume  larger  than  that  of  the  convex  body.   If  the  differ- 
ences are  analyzed,  it  can  be  shown  that  the  form 


Ip  3  I     ..  3P6 


+  cp  (A4-6) 


9P  'CB     3P  'EHS 


yields  excellent  agreement,  where  the  superscripts  are 
CB  for  convex  body  and  EHS  for  equivalent  hard  sphere. 
This  may  explain  much  of  the  success  that  Mathias  (1979) 
had  in  his  correlation  of  dense  fluid  properties. 

The  question  of  whether  the  RISM  approximation  employed 
in  this  work  can  yield  good  agreement  for  hard  body  systems 
is  more  difficult  to  answer  due  to  lack  of  data.   The  direct 
correlation  values  given  by  Lowden  and  Chandler  (1973) 
for  a  hard  diatomic  system  have  been  integrated  to  yield 
the  thermodynamic  derivations.   These  are  shown  in  Figure 
A4-3.   Chen  and  Steele  (1971)  have  solved  the  Ornstein- 
Zernike  equation  for  this  same  system  using  a  Percus- 
Yevick  approximation. 
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Figure  A4-2.   Dimensionless  compressibility  versus 
dimensionless  density  for  hard  core 
systems . 
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Figure  A4-3.   Dimensionless  compressibility  versus  density 
for  homonuclear  diatomics  calculated  using 
RISM  direct  correlation  functions. 
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A  comparison  of  these  results  is  shown  in  Figure  A4-4. 
The  RISM  theory  offers  a  good  approximation.   It  should 
be  noted  that  Chen  and  Steele's  results  came  from  the  com- 
pressibility relation  using  the  P-Y  approximation.   Thus, 
they  would  be  expected  to  overestimate  the  exact  result 
for  a  diatomic.   This  then  offers  even  greater  hope  for 
the  RISM  theory.   Figure  4-2  in  Chapter  4  also  shows  a 
comparison  of  the  RISM  result  for  this  system  with  results 
from  Monte  Carlo  calculations;  again  the  comparison  is 
favorable. 

These  results  suggest  that  the  RISM  formulation  should 
be  more  accurate  for  the  representation  of  the  hard  core 
system  than  use  of  a  spherical  reference,  but  that  if  prop- 
erly formulated,  the  spherical  systems  can  also  work  well. 

Effects  of  Long  Range  Interactions 

For  polar  molecules  the  effects  of  the  multipole  moments 
on  thermodynamic  properties  can  be  large  in  the  vapor  phase. 
However,  in  the  liquid  phase  the  harsh  repulsive  forces 
seem  to  dominate.   A  prime  example  of  this  behavior  is 
shown  in  Figure  A4-5  for  argon  and  water.   The  reducing 
parameters  for  water  were  chosen  to  yield  the  optimum  cor- 
respondence for  these  two  substances.   It  can  be  seen  that 
these  two  different  fluids  behave  similarly  in  the  dense 
fluid  region  while  being  highly  dissimilar  in  the  vapor 
region. 
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Figure  A4-4.   Comparison  of  RISM  and  integral  equation 
results  for  homonuclear  diatomics  with 
L/o  =  0.6. 
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Figure  A4-5. 


Comparison  of  the  Bulk  Modulus  for  Argon 
(lines)  and  water  (o).   e  and  o  for  argon 
are  from  virial  coefficient  calculations 
and  e  and  o  for  water  were  determined  for 
best  correspondence. 
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Gubbins  and  O'Connell  (1974)  have  analyzed   this 
behavior  for  multipolar  potentials  and  spherical  molecules. 
They  find 


3P[ 
3p 


3PB 
3p 


[1  -  Gx/[3PB/3p 


]] 


!A4-7) 


where. the  subscript  o  represents  the  properties  of  a 
spherical  reference  system.   Using  a  Lennard-Jones  reference 
system  the  quantity  G   was  found  to  be  approximately  -1 
for  a  dipolar  interaction.   In  the  dense  fluid  state 


0.05  so  that  the  anisotropic  forces  contribute  approxi- 


3PS 

3p 

mately  5%  to  the  value  of  3P3/3p.   This  suggests  that  a 

low  order  perturbation  expansion,  as  used  in  this  work, 
should  be  adequate. 

Summary 

This  appendix  has  been  dedicated  to  a  study  of 
nonsphericity  effects  on  thermodynamic  properties.   Two 
separate  considerations  were  examined--the  shape  of  the 
molecule  and  anisotropic  long-range  forces.   It  was  shown 
that  the  nonspherical  shape  could  be  adequately  represented 
using  the  RISM  theory  or  even  with  the  properly  chosen 
spherical  reference.   The  anisotropic  forces  can  be  a  notice- 
able effect  but  the  low  order  perturbations  should  be 
adequate  for  the  systems  of  interest. 


APPENDIX  5 
HARD  SPHERE  PROPERTIES 


Exact  expressions  for  the  correlation  functions  and 
thermodynamic  properties  of  hard  sphere  mixtures  are  not 
available.   However,  good  approximate  results  can  be  obtained 
using  the  Percus-Yevick  aproximation.   Certain  aspects 
of  this  approximation  and  its  ramifications  on  the  form 
of  the  hard  sphere  equation  of  state  are  developed  below. 

To  determine  the  properties  of  any  system  of  molecules, 
one  begins  the  analysis  with  the  Ornstein-Zernike  (1914) 
equation 


h.  .(r)  =  c  .(r)  +  I      p.  J  ds  c,  (r-s)  h.  .  ( s )       (A5-1) 

1]  1]         ,      K  IK  K  J 


The  form  of  the  Ornstein-Zernike  equation  is  valid  for 
any  system  of  molecules,  acting  with  spherically  symmetric 
intermolecular  potentials,  not  in  an  external  field.   This 
represents  a  set  of  N(N+l)/2  coupled  integral  equations 
that  contain  N(N+1)  unknown  functions,  for  an  N  component 
system.   For  a  mixture  of  hard  spheres  there  are  a  compli- 
mentary set  of  exact  conditions 


h.  .  (r)  =  -1  for  /r/  <  o .  .  (A5-2 
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Ill 

which  represent  the  impenetrability  of  the  hard  cores. 
For  the  cases  of  interest  the  characteristic  length  o 
is  given  by 


o.  .  =  -  (a.  +  o  .  )  (A5-3) 

ID    2    i     D 


where  o.  is  the  diameter  of  molecule  i.   Percus  and  Yevick 

l 

(1958)  have  proposed  a  set  of  approximate  closure  relations 
for  the  direct  correlation  functions  that  when  combined 
with  equations  A5-1  and  A5-2  yield  a  set  of  equations  that 
can  be  solved  for  the  direct  correlation  functions  and 
the  equation  of  state.   These  approximate  relations  for 
a  hard  sphere  system  are 


c.  .  (r)  -  0       for  /r/  >  a.  .  (A5-4) 


Werthein  (1963)  and  Thiele  (1963)  offered,  indepen- 
dently, solutions  to  this  problem  for  pure  components. 
Later,  Lebowitz  (1964)  gave  the  solution  for  a  binary  mix- 
ture.  Baxter  (1970)  has  shown  an  elegant  approach  for 
the  solution  for  multicomponent  systems. 

Once  the  correlation  functions  are  known  for  the  hard 
sphere,  the  equation  of  state  can  be  determined  in  two 
obvious  ways.   The  so-called  compressibility  result  is 
found  from 
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|££)    =1-1  X.P,  I  dr  c   (?)  (A5-5) 

9p   T,X        i   X  D        ±J 


with  use  made  of  the  ideal  gas  limit  to  express  this  as 
a  complete  equation  of  state.   The  other  form  of  equation 
of  state,  known  as  the  virial  results,  is  found  from  the 
virial  theorem 


PB  «  P  -  (f)  *  11    PiPj/0  dr  r   ("d^Kj  (r)       (A5_6) 


Another,  more  roundabout,  approach  to  obtain  an  equation 
of  state  is  possible.  One  may  first  obtain  the  internal 
energy  and  then  use  thermodynamic  identities  to  find  the 
pressure . 

Because  the  correlation  functions  determined  with 
use  of  the  Percus-Yevick  approximation  are  not  exact,  the 
different  routes  to  the  equation  of  state  all  yield  different 
results.   For  pure  components  the  closed  form  results  avail- 
able are 


2        3 
Compressibility:   Z  =  (1+n+n  )/(l-n)  (A5-7a 


Virial:  Z  =  ( 1+ n+ n2-3 n3 ) /( 1- n) 3     (A5-7b) 


where  n,  the  packing  fraction,  is  found  as 


113 


n  =  ifhpo3 


These  results  are  known  to  be  upper  (compressibility)  and 
lower  (virial)  bounds  to  the  exact  hard  sphere  results 
as  determined  from  computer  simultation  (Reed  and  Gubbins, 
1973). 

By  comparison  of  the  density  series  expansions  of 
the  exact  and  approximate  results  for  pure  hard  spheres 
Carnahan  and  Starling  (1969)  developed  an  expression  which 
very  accurately  represents  the  hard  sphere  equation  of 
state 


z  .  (i+n+n2-n3)/(l-n)3  (A5-8) 


A  generalization  of  all  of  these  results  can  be  written 
as 


z  =  (1+n+n2+an3)/(l-n)3  (A5-9) 


where  different  values  of  the  parameters  a  can  yield  any 
of  the  aforementioned  forms. 

There  are  several  ways  that  a  value  of  a  could  be 
chosen.   One  reasonable  approach  would  be  to  use  a  to  match 
the  known  virial  coefficients  for  hard  sphere  fluids.  Table 
A5-1  shows  the  values  of  a  required  to  reproduce  the  known 
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TABLE  A5-1 

VALUES  OF  ota  (EQUATION  A5-9)  REQUIRED  TO  MATCH  KNOWN*3 
VIRIAL  COEFFICIENTS  FOR  HARD  SPHERES 


Virial  Coefficient  a 

4th  -0.6384 

5th  -0.9211  +  0.0256 

6th  -1. 0789  ±  0.0683 

7th  -0.7475  ±  0.1638 

The  second  and  third  virial  coefficients  are 
independent  of  ex. 

The  fourth  virial  coefficient  is  evaluated  exactly, 
whereas  the  fifth  and  higher  virial  coefficients  have  been 
determined  using  Monte  Carlo  integration,  and  the  errors 
indicated  are  based  on  the  statistical  error  of  those 
calculations . 
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virial  coefficients  (Hansen  and  McDonald,  1976)  for  hard 
spheres.   These  results  show  that  the  Carnahan-Starling 
equation  (A5-8 ) ,  for  which  a   =  -1,  gives  a  reasonable 
approximation  for  the  known  virial  coefficients.   In  this 
work  a  will  be  treated  as  a  parameter  determined  through 
analysis  of  experimental  data  as  described  in  Chapter  5. 
An  important  consideration  in  the  choice  of  a  value 
for  a  is  that  the  hard  sphere  systems  show  no  vapor-liquid 
phase  transition.   This  implies  that 


2* >  o  o  <  n  <  1  (A5-10) 

dn 


Using  equation  (A5-9)  it  is  possible  to  predict  a 
maximum  in  Z  for  certain  values  of  a.   This  is  actually 

the  case  for  all  negative  values  of  a.   The  location  of 

HS 
the  maximum  in  Z    is  given  by 


±  /  62-4l 


max  _  -P  ±  j     P"-4P  (A5-11 


;ith  8  determined  as 

6  =  4/(3a  +  1)  (A5-12) 

Analysis  of  these  results  shows  that  any  value  of  a 

LJC 

less  than  -3  will  yield  a  maximum  in  Z    for  values  of  n 
less  than  1.   However,  on  a  purely  physical  basis  n 
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should  never  exceed  that  for  closest  packing  of  spheres 
which  is  approximately  0.74.  For  the  value  of  a  used  in 
this  work,  -4.2,  equatin  A5-9  will  show  a  maximum  in  the 
compressibility  factor  at  n  equal  to  0.784.  Thus,  these 
results  should  probably  not  be  used  for  n  values  greater 
than  0.75. 

To  extend  these  results  to  mixtures  a    can  be  treated 
as  a  constant  that  acts  as  an  interpolating  parameter  between 
the  compressibility  theorem  (denoted  as  PYC )  and  virial 
theorem  (denoted  as  PYV)  results.   The,  compressibility 
factors  and  chemical  potentials  are  found  using 

ZHS  =  ZPYC  +  (0.)  [ZPYC  _  ZPYV]  (A5-13a) 


5yHS  =  ByPYC  +  (SL)  [ByPYC  _  ByPYC  }  (A5-13b 

i       i       J      i        J- 


The  complete  expressions  for  all  of  the  required  quantities 


ZPYC  =  ±    [_fo  +  JJ1J2   +  3g2 ]  (A5-14a) 

'o       1_C3  (1-C3)2  (I-C3)3 

ZPYC    _    ZPYV  =       3^2  (A5-14b) 

^o(1^3)3 

6yfC    =    BUi    "    ^n    (I-C3)    +   o3    ,Q    ZPYC 


2       9q  2   2 
:  _._-_^_  +    1  ^   +  _li L_ (A5-14c 


3°i^2   ,  3°x  J!   ,  9°i  g 
(!-£.,)    (!-£,)     2(1-5 


2 
3/     u-s3-      «.»*  s3) 
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2 
,yPYC    _    ByPYVj    z    _±_2    ,n(1-C3)[3C22o.     -    1/2] 


K. 


30. 3    ,    3 

+  — i ^    [    6    +  -     -, 

2(l-?3)3  C3 


2(1-5)         4(1-    | 


1-C3) 


C3 


where 


3o.z£ 

+       X      2   2 
2(l-53) 


'"* 


(A5-14d) 


?£    =    6    ?    Pi°i 


(A5-14e 


APPENDIX  6 
COMPUTER  PROGRAMS 


AJGJST  26,   1S31         1  1".  29  NEROC   —    SYSTEM  SUPPORT  UTILITIES 


C  PROGRAM  HARD  SPHER.- 

C 

C  ThI5     IS     DESIGNED    TO     BE     A     MULTIPURPOSE    PROGRAM     FOR 

C  FITTING     DATA     TO    THE     PERTURBED    HARD     SPHERE    FORM 

C  PROPOSED     FCR     ThE    DIRECT     CORRELATION    FUNCTION     INTEGRALS 

C  THERE     ARE     T*C     GENERAL     OPTIONS     AVAILABLE    TO    THE    USER 

C  NOPT=  1 HERE     THE     INPUT    DATA     IS     THE     TYPE    CF    HARD     SPHERE 

C  ECUATICN     TO    aE    USED    ANO     THEN     ISOTHERMAL     SETS 

C  CF     VOLUME     AND     DCF I     DATA     ANO    A     VALUE    OF    THE     hARD 

C  SPHERE     01 AMETER(DIMENSIONLESS)     FOR     THE     ISCTHERM 

C  THE     FRCGRAX    WILL     CALCULATE    THE    HARD    SPHERE    OCFI 

C  AND     THE    CORRESPONDING     VIRIAL     COEFFICIENT 

C  NOFT=2 WITH    THIS     OPTION     AN    OPTIMAL    VALUE     OF     THE 

•C  HARD     SPHERE    DIAMETER    WILL     BE     FOUND    FCR    AN 

C  ISOTHERM     BY    MINIMIZING    THE     SUM    SQUARED 

C  DEVIATIONS     OF    THE     VIRIAL     COEFFICIENTS    FRCM 

C  THE     MEAN    ALCNG    THE     IaOTHERM 

C  AGAIN     SEVERAL     FORMS    OF    HARD     SPHERE     EQUATION 

C  OF     STATE     ARE    AVAILABLE 

C 

C  ThE     INPUT     CATA     IS     SLIGHTLY     DIFFERENT    FOR    THE    DIFFERENT 

Z  JPTIONS.         IN     eCTH     CASES    ThE    FIRST     CARD     SHOULD     CONTAIN    THE 

Z  COMPONENT     NAME     AND     CRITICAL    PROPERTIES.        THE    NEXT     CARD 

C  CCNTAINS    Th=    VALUE     CF    THE     PARAMETER    AL°HA    THAT     DETERMINES 

C  THE     TY^E     OF     HARD     SPHERE     EQUATION    USED. 
C  AL^HA  EQUATION 

C  0  PY  C 

C  -3  PYV 

C  -t  CS 

C  -3TC-1  SOFTER     THAN     CS 

C  -  1     TO     0  HARDER     THAN    CS 

C  THE     NEXT     CARC     CONTAINS     NOPT     =1     FOR     CALCULATION    AND 

C  -2    C0P     REGRESSION 

C  NEXT     CARD     IS    THE    FIRST     TEMPERATURE 

C  Z.ST2R     0     IF     THE     END     OF     AN     ISOTHERM 

C  ENTHH     NUMBER     LARGER     THAN     10**4    FOR     A    NEW    COMPONENT 

C  THEN     SETS     CF     V     AND     C    EIGHT     TO     A     CARD 

C  FINALLY      IF      NCFT     =1      ENTER     XMAX.XMIN     AND     XINT 

C  IF     J0PT=2    ENTER     INITIAL     GUESS    FOR     X 
C 

C  ALL    CALCULATIONS     IN    DOUBLE     PRECISION 
IM°LICIT     REAL»£     (A-H, O-Zj 

CI  KENS  ION     Vl(e).CIO),  FOL(  10  ).  IN  A  {  3  )  .  V  (  30  )  ,  C{  30),  DEL  (  30). 
»       E( 30) ,NA«E(5>  ,CHS( 30 )  ,A(  7) 
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C  RESTRICTED     TC     30     POINTS     ON     AN     ISOTHERM 

CJMvlOh     /TITLE/    NAME.TC.VC 

CCl»MON    /WRITER/    NOPT 
C  ENTER     N4ME     AND     CRITICAL     PROPERTIES 

1  30  READ ( 5.  1 . END=9  999JNAME.TC, VC 

1  FORMAT <5A4,2F10.0) 

C  JRITE     HEADER     AND     THEN     ECHO     NAME     AND    TC     ANC    VC 

*  SITE  (6,  2) 

2  PufiM4T(  I  OX  ,•  ANALYSI  S    OF     DIRECT     CORRELATION    FUNCTION' 
*.  •     M3DEL     USING    HARD     SPHERE     REFERENCE    SYSTEM*,///) 

«RITE(6,3)NAVE,TC.VC 

2  F  JRV|AT(  1  OX, 'DATA     FOR     •  ,5A4  ./// . 1  OX ,  • C RI TI CAL    PROPERTIES" 

*  /.  1  OX.  'CRITICAL    TEMPERATURE    =  • , FS. 3.  • K'  . / , 1  OX, 

*  'CRITICAL     VCLUME     =' , F9.4 . • CC/MOLE ' ./// ) 

C  ENT£R     OPTICN    FOR    CALCULATION    OR     REGRESSION 

REA0(5,6  )NCPT 
£  FORMA T( 12) 

C  ECHO     TY°E    OF    CALUCLATION    PERFORMED 

GC     TO     ( 1 00C  ,2000) • NCPT 
1000        *RITE(6.7) 
7  FC«MAT< I  OX,  'PROGRAM     IN     CALCUL AT IONAL    MODE', ///I    . 

3C    To    3000 
2000        »RITE(£.£) 

3  FCPM AT {1  OX, • PROGRAM  IN  REGRESSION  MODE',///) 

C       SET  COLNTER  FCR  NUM3ER  CF  CATA  POINTS  CN  AN  ISOTHERM 
3T0J   NP=0 

C       TOL")  USED  TC  DECIDE  IF  A  NEW  ISOTHERM  IS  ENCOUNTERED 
TJL )=0. 030 
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;JGJ5T     26.      1921  li:29  NERDC       SYSTEM     SUPPORT     UTILITIES        


C  ENTER     THE     TEMPERATURE 

210  REAOl  5.  4)T 

i  FCRMAT (Fl 0 .0  ) 

C  CHECK     TO     SEE     IF     THE     END     GF     AN     SCTHERM 

IF<T.LT.  1.000 )     GO     TO     400 
C  CHECK     TC     SEE     IF    END    OF     A     CCMPONENT     DATA    SET 

IF( T.GT.  1. CDC4)     50     TO     100 
C  INPUT     VOLUME     ANO    OCFI 

RcAO (5t9)VI  .CI 

9  FQRMAT( 8F1C.0) 

C  CHECK     TO     SEE      IF    A    NEW     ISOTHERM 

Ir(NO.EQ.O)     GO     TO     300 

IFOABS(  T-TOLO  )  .GT.  0.5D0)     NP=0 
C  N3»     PACK     CATA     FCfi     THE     ISOTHERM     ANC     CALCULATE    NP 

300  00     10     1=  If  £ 

IF(V  I  (  I  >  .LT.l  .000  )     GO     TO      10 

NP=NP+1 

V ( S3 ) =  V  I  (  I  ) 

C (NP)  =CI  (  I  1 

10  CONTINUE 
TGLO=T 

C       kcTJRN  FOR  NEXT  ISOTHERM 

go  t:;  200 
c     decide  what  to  read  next 

400  GC     TO     (ll*12)tKCPT 

11  S£AD( 5, 1 01 ) XMAX.XM IN.XI NT 

1  11  FORMAT (3F10.0  ) 

R5AD( 5. 1 05 )NA, A 
105     FOKMAT< I  2.  EX.  7F1 C. 0) 

CC  20  K= I . NA 

AL=A( K) 

*RITE(b. 5)AL 
5       FOSMA T(  1  OX,  • VALUE  OF  ALPHA  IS  * • FB . 2. // • lOXt 

*  UL3HA  =  C   IS  PYC  . /.I  OX,  •  ALFHA  =  -3   IS  PYV'./.IOX, 

*  «ALPHA=-1   IS  CS1,///) 

C  CALCJLATE     NLMBER     OF     VALUES     CF     X     TC     eE    TRIED 

NX=l      ♦     < XMAX-XMIN) /XIN T 
C  EEGIN     CALCULATIONS     F3R     DIFFERENT     X     VALUES 

DC    2  0    J=l  .NX 

<=     XM IN+XINT*DFLOAT{ J- 1) 

CALL     FUN(CfV.AL.X,F,VC,NP, CHS, E.OEL) 

CALL     OUT(XtC.V,F. TOLO, NP.CH S» E , DEL) 

2  0  CONTINUE 

C  RETURN    FCR     NE  »     ISOTHERM 
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3C    TO     20  0 
C  READ     CATA     FOR     REGRESSION 

12  REAQ(5.  A  )X 

REAO(  5.  1  05  )  h>A.  A 

3  3     50     K=  It  NA 

AL  =  A (K) 

•  RITE (C  .5) A  I 

WR  ITE(6.   126  )X 
126  FCFMAT  (1  OX  .  •  IN  IT  IAL     GUESS     OF    X    =  •  i  F  12.  5.///  ) 

Z  3EGIN     REGRESSION    USING     RQUADO 

1TEST=2 
?100        CALL    BCUADCd  TEST,  X  ,F,  100.  1  .OC-6.  1  .00-6,0  .1D0.  RDLi  INA  ) 

GO    TO     (  1  IOC,  1200,  1 300)  • I  TEST 
C  EVALUATION     OF     SSC 

110J        CALL    FUN(C.  V.AL.X.F  ,  VC  ,  NP . CHS » E . DEL  ) 

Gu     TO     2100 

C  REGRESSION     COMPLETE PRINT     RESULTS 

120J        CALL     OLT( X.C. V.F.TOLO.NP .CHS.E.OEL) 
C  hicTJRN    FOR     NEW     ISOTHERM 

GC    TO    50 
C  ERROR     IN     REGRESSION 

1300         CALL     EHR?P (X.F, ITEST.TCLD) 
C  RETURN    FOR     NE*»     TEMPERATURE 

50  CONTINUE 

GC  TO  200 
SS!»i   STOP 

ENC 

SuEROLTINE  CALCI AL . X.V.E.CHS.VC ) 
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THIS     ROUTINE    CALCULATES    THE    DCFI     FOR    THE    GENERAL    HARD 
SPhcRE     EQUATION    CF     STATE 

ALL     CALCULATIONS     IN     DOUBLE     FRECISICN 
IMPLICIT     REAL*3     (A-H.O-Z) 

CALCJLATE     THE     PACKING     FRACTION     ANO     CCF I 
£=X*VC/V 

CHS=  (  (I  . ODO+AL  )*E**4-4.0D0*(  I  •  OD  0*  AD  *E*»  3*2.  0O0»E*E- 
»        8.  0D0«E) /<1.  0C)0-E)**4 
RETURN 
ENC 

SuBMJLTI NE    FwN(C.V.AL.XtFtVC.NPtChS.EtOEL) 
THIS     SUBROUTINE     CALLS     CALC     TO    EVALLATE     THE    HSDCFI     FOR 
THE     VALUE    CF     X     SENT     FRCM    RCUADC     ANO    THEN    FINDS     THE 
VIRIAL     TER*«    FOR    EACH     OENSITY     AND     AVERAGES.        THE 
OUTPUT     IS     THE     SUM     SQUARED     DEVIATIONS    FROM     THIS    MEAN 
ALL     CALCULATIONS     IN     DOUBLE     PRECISION 
IMPLICIT     REAL*6     (A-H.O-Z) 

DIKE NS ION     C(1).V(1).CHS(1).E< I). DEL (I) 
CCfMON     /ASiG/DMEAN 
INITIALIZE     MEAN     TO     ZERO 
CE=O.0CO 
DO     13     1=1 iNP 

CALL     CALC(  *L.X,V(  I  ).E(   I  )  .  CHS(  I  )  .  VC  J 
JEL( I  ) =(  C<  I  >-CHS(I)  )*V  (I  ) 
J3=D3+DEL(  I  ) 
FINO     AVERACE 
D*EAN=DB/NP 
INITIALIZE     SSD     TC     0 
F  =  0.OC0 
DO     2  0     I=1.NP 
Fl=(DEL(   I)-DMEAN)**2 
FsF  +  F 1 /V( I  ) 
RETURN 
ENC 

SLdROUTINE     CLT     ( X , C  .V , F ,T . NP . CHS .£. DEL ) 
THIS     SUBROUTINE     IS     USED     TO     PROVIDE     OUTPUT    FROM    THE 
REGRESSION     MODE    OF     HARD     SPHERE 
NO     CALCULATIONS    ARE    PERFORMED 
IMPLICIT     REAL*3     (A-H.O-Z) 

DI  KENS  I  ON     CO  )  tV  (1  I  tCHS  (  1  )  .  E  ( 1  )  .  DEL  ( 1  )•  NAME  (5  ) 
COMMON     /AVG/DMEAN 
COMMON     /TITLE/    NAME.TC.VC 
COMMON     /WRITER/    NOPT 
CALCULATE     AND     OUTPUT     THE    RECUCED     TEMPERATURE 
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*  FCPMAT(10X,Fc.3.2X,F7.4,2X.F7.4.2XtFe.2t2X.F8. 3t2X,FlC.A. 

*  2X.F3.  3.2X.F15.4  > 
3=S/NP 

WRITE (6,5  >  S 
?  FCR^ATt  10X,  'AVERAGE     ABSOLUTE    DEVIATION     IN    C     IS     '^10.4,///) 

RETURN 

£NO 

SUaROUTINE    EPROR(X ,F, I ,T) 
:  THIS     SUBROUTINE     WILL    PFINT     ERROR     MESSAGES    RESULTING 

C  FRO*     PROBLEMS     ENCOUNTERED     IN    RQUAOD 

IMPLICIT     REAL*3     CA-h,0-Z) 

I  =  1-2 

GO     TO     (  1  .  2  )  .  I 
1  *RITEC6. 3JT.X.F 

3  FQRMUI  10X,  'ERROR     IN    RQUADD    DUE    TC     RCJNCOFF"  . /// , 1  OX, 

»        "TEMPERATURE    =     • . F  12 . 4, //, 1 OXt  • LA  ST     VALUE    OF     X     =     '.F12.4, 

*  //.lOX.'FUNCTI CN     =     »,D15.5,///) 
3c  1  JRN 

1  «3ITE{6. 4)T,X.F 

4  FOH*AT< 10X. 'ERROR TOO    MANY     ITERATIONS    AT     T= • ,F 1 2 .5. // t 

*  10X,  'VALUE     CF     X     I S« , F12.4.2X. • *  I TH    F      =•  ,D 1 5. 5. /// ) 
RETURN 

END 

3NU     JF    FIL-E     ON     SYSIN. 
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Tfi-T/TC 

iKH  ITE  16.D1.TR 
1  FCfi MAT ( I  OX, • TEMPERATURE    = • , F9 .3 • 5X , «T/TC    ='.F10.4, ///) 

C  OUTPUT    HAD     SPHERE     D  I AMETER . AVERAGE     VIRIAL    COEFFICIENT    AND 

C  5UM    SCUAREC     CEVIATI3N 

3C    TO     (  ICC. 200) .NCFT 
>00  *RITE( 6. 2) X.CMEAN, F 

'  FGR»4AT  <  1  OX,  'THE    OPTIMUM     HARD    SPHERE    DIAMETER    WAS    FOUNO     TO     HE' 

*  Fl 2.5, 3X,  'THIS    YIELDS    AN    AVERAGE     VIRIAL     COEFFICIENT     OF', 

»       Ft  2  .5,//.  10X,  'THE     S J M     SQUARED     DEVIATION    FROM     THE     MEAN     =•. 

*  315.5.///) 
30     TO     30  0 

100  4R ITE<6, 29  )X, CMEAN. F 

2?  FuRMAK  1  OX,  •  VALUE     CF    X      =', 

*  F  12.5.  2X,  'THIS     YIELDS    AN     AVERAGE     VIRIAL    COEFFICIENT     OF', 

*  F12.5.  //  .10X,»  THE     SUM    SQUARED    DEVIATION    FROM     THE     MEAN    ='. 

*  Dl 5.5. ///) 

C  *RITE     HEADERS 

300  «RITE(6.3) 

3  FOR-4AT(  IOX,  'VOLUME  RED     DEN  ETA  DCFI  HSDCFI', 

*  '  VIRI/L  CCALC  ERR     IN    C  •.////) 
S=C. ODC 

DO     10     1=  I,  NP 
RC  =  VC/V (  I  ) 

CCAL=CHS(I) +OMEAN/ V<  I) 
PE=C< I )-CCAL 
S=S+S»eS (PE) 
10  «RITE(e,4)  V(  I )  .RD.E  (  I)  .C(  I  )  .CHS  (I  )  .DEL (I)  .CCAL.PE 
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//OLI  ST     JOB     (I  00  6,064  3.5.  9. .2061.  2. N.  48).  'TELOTTE  •  •  CLASS=  A,  TYPRLN=COPlr        JC 

/*KJUTE     ORINT     LOCAL 

//    EXEC    F3RTG-C 

//FuRT.SYSIN        CC        * 

C  PROGRAM     TO     CALCULATE    P-V-T     PRCFERTIES     CF    ARGCN 

C  ThIS     PROGRAM    USES     ThE     GENERALIZED     STRCtiRIDGE    EQUATION 

C  DEVELCPED    BY     TfcU.LEE    AND    STARLING    FCR     CALCULATIONS 

C  FOR    LIQUID     ARGON 

C  EQUATION    CF     STATE     CALCULATIONS     USING    THE     STROBRIDGE    ECUATICN 

C  CALCULATIONS    DCNE     *  I TH     01 MENSI ONLESS    VARIABLES 

C  Tt-E     LOCATICN     IN     THE     PROGRAM    TO    BE     CHANGED    FOR 

C  LSE     JF     DIFFERENT    RECUCING     PARAMETERS    WILL     BE     INDICATED 

C  CONSTANTS    FOR     THE    DIFFERENT     SUBSTANCES     ARE    TC     BE    ENTERED 

C  IN     A     CATA     STATEVENT 

IMPLICIT     REAL*b     (A-H.C-2) 

0 IMENS ION    A(  16 ) 

CATA  A/1  .3  10  24.-3.8  063  6. -2.  3723  6. -C. 798 672.0. 198761. 

*  1.47  014.-C.766367.2.  15465.5.7  54  29  .6.7  822.-9.94904. 

*  -1 5.6  162, £6.643.  1  £ . S 27, 9 . C47 55 . 6. 63232/ 
READ (5.  I >NA*E1 .NAME2 

1  F0RMATC2A4) 

*RIT£(6.2>NAME1.NAME2 
Z  FQRnIAT<  1  OX  ,•  ECUATICN    OF    STATE     CALCULATIONS     FOR     «.2A4.////) 

READ(  5.  3 )TMIN,TMAX.  TINT 

REAa(5,3)RMN.PWAX.  FINT 

3  FCRMATI3F1 C. C) 
JCCJNT=(TMfiX-TMIN)/TINT+l 
KCCJNT=(RV  AX-RMIN)/R  INT+1 
DO      10     I=1.JCCLNT 
T=TMIN+OFLOAT ( 1-1) *TIMT 

Z  HERE     IT    HAS     BEEN    ASSUMED    THAT     The     kAX. MIN     AND     INTERVALS     FOR 

C  THE     TEMPERATURE     ITERATIONS    HAS    EEEN    ENTERED     IN    OI MENS  I CNLESS 

C  ORIAELE3 

C  CHANGE    THE     NEXT     STATEMENT     TO    CHANCE    THE    FORM    UF    DI MENS ICNL ESS 

C  TEMPERATURE.        WE     HAVE     USED     HERE     INPLT     AS    T/TC     AND     THE 

C  CALCULATICNS     »ILL     BE     DCNE     *  IT  h     K*T/EPSILON 

C  THE    RELATICNSHIP     IS    EP SI LON/K=TC/ 1 . 25 33 

TS=l ,2593CC*T 

iRI  TE  (6  ,4)  T  ,TS 

4  FJR4AT(  10X.   'REDUCED     TE MPER ATURE=  •  .F 5.4 . 5X. •  T ST AR  =  «  . F9 .4  ./ //  . 

*  lOX.'RECUCEC     CENS  ITY «,5X,  •  R HO     STAR'.IOX,'  Z  '.ex. 

*  ■         1-C        '  .///) 

C  NCW     DEFINE    TEMPERATURE     DEPENDENT     FLNCTICNS    REQUIRED 

TI=1 . ODO/TS 

F1=A(1)+TI*(A(2)+TI*(A(3)+TI*(A(4)+A(5)*TI*TI))) 
F2=A(6  )*A(  7  )*T  I 
F3^=A(  6) 
F4=TI  *TI*TI*(A(9)  +  Tl*(A(10)+A(ll)*TI)) 


127 


F5=TI*TI*TI*(A(12)+TI*(A 
F6  =  A(   IS ) *TI 

CC     11     J=li  KCCUNT 

R  =  RMI N+DFLCATI J- 1 ) * R I N T 

RS=3.2  189CC+R 
C  TmE     LAST     STATEMENT     SHOUL 

C  DENSITY     VARIABLES.         THE 

C  Tl-E     PROGRAM    USES    RS=*H3* 

C  5lGrtA**3=0. 3189/RHCC 

C  »ILL     NOW     DEFINE    DENSITY 

R2  =  RS*RS 

S3-RS»RS*fiJ 

r*4=3S*RS*RS  »RS 

P5=RS*RS*RS*FS*RS 

0  l=-A  (  16)»R2 

El=OEXP( CI  ) 
C  NUl     CALCULATE     TME     Z     ANC 

C  *ITH    RES°ECT    TC     THE     DENS 

F=l  .0  CO  +F1  «fiS  +F2*R2+F3*R 

3F=F 1+2. 0DC*F2*RS+3.0DO* 

♦  -2.0C0*A(  16  >*R3)*FS*E1 

•  +  3. 0C0*F6*R4 

qmc=r  s*df+f 

11  «rtITE(6. 5)R.RS.F,0MC 

5  FCfiMA  T (  1  3X  ,F5  .4  ,6X  ,  F9  .  4  . 

10  *RlTE(6t6) 

5  FCRMAT (1HO ) 

SIC? 

END 


(13)+* (14  )*TI)  ) 


C     EE    CJ-ANCE3    TO    USE    DIFFERENT 
INPLT     VALUES     WERE     R=RHC/RHOC 
SIGMA**3     *ITH     THE    CORRELATION 

DEPENDENT    FUNCTIONS 


DERIVATIVE    CF    The    Z 

ITY 

3+F4*R2»El+F5*R4*El+F6*R5 

F3*R2+F4*E1 *(2 .ODO*RS 

*(  4.0DO*R2-2.0DC*A(  16)*R5) 


10X.F9.4.  12X.F9.4  ) 


Jc32     JGB    STATISTICS 


2b    ALG     £1     JOB     EXECUTION     DATE 


73     CARDS     CEAC 


0     SYSOL1     PRINT     RECORDS 
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C  PUPPOSE:         TC     GENERATE     P-V-T     PROPERTIES     OF    LIQUID    METHANE 

C  THIS     PROGRAM     LSES    THE     EQUATION     CF     STATE    OF     MCLLERUP 

IMPLICIT     REAL*8     (A-h.O-Z) 

3SNERIC 

COMMON     A1,A£,A3.A4,B1,B2.B3.8*»C1»C2.C3.C4.C5.C6,D1.D.2.D3.C4,05 

COMMON/ONE^AC-.  AL.BB.BE.DE.  GK,     TTRP.PTRP.OTRP.     TCRT t PCRT.DCRT 

CCfrM0N/THRE/3PDT.D2PDT2.DLPDLR.0TS0S.DTHOR.DPSIDR,     XB  1  . XB 2 . XO 1 . XD2 

COMMON/FIVE/     D,  T.  DEN.DPDDi     E.H.S.     CV.CP.Ii.     WK 

CCI»MGN/E  IGt-T/     ODLDT 

COMMON/FOLF/FPT.DFDX.DFDRl 

CINENSION     PP( 16) i VV(  16 ) tOMCC  16) 

CALL     NETHAh 
100  REA:>(  5.  l.ENO=S9)  T 

1  FOR  *AT  (FIO  .0 ) 
*R  I  Tc  (  6  t  2  I 

2  FORMAT(  lHlt  10X.  'METHANE    PROPERTIES     FROM    MOLLERUP* , •// ) 
TP=T/TCRT 

aR IT:( 6.3) T  ,TR 

3  FGS«1ATl  10X,  'TEMPERATURE     = • . F 1 0. 3. •/ , 

i        •  REDUCED    TEMPERATURE    =•  .Fl 0  .*./// ) 

MR  ITE  (  6.  4) 

4  FCP-MAT  <  1  OX.  •  PRESSURE  VOLUME  1-C«.///) 
BAS~  = 15. 22EDC 

UINC=  1  .0  15CC 
CC     10     1=  1  .  16 

VV(I)     =     3AJE     +     DFLOAT( 1-1 ) *D1 NC 
PP(  I  )  =  0°DRF (T. VVt I  )  ) 
UMCI  I  >=OLPQLP*PPU  > /VV  <I)/T/GK 
VV( I ) = 1. OD c/VV( I ) 
Pt;U)=PP(I)/l.O132  5D0 
10  mrtlTE  (  6.  5)  FP(  I  >  .  VV  (I  )  .  C  MC  (  I  ) 

5  F0RMAT(  1  OX, 3F 15.4) 
»PI TE  (7 t6)T 

»R IT^( 7. 6)PP. VV.GMC 
*K  IT=(7.  6  )T 

6  FCfiMATOFl  0.3) 
GG     tj      100 

35  STCP 

ENO 

SUBROUTINE     VETHAN 
C  NCTc     THAT     PRESSURES     ARE     IN     BARS.      1     ATM     =     1.01325    BAR. 

C  ONE    6AR-LITEP-     =     100      JOULES. 

C* 

IMPLICIT     PEAL*S     (A-h.C-Z) 
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GENcR  1C 

CCM^ON    Al,A2,A3iA4,31iE2.SJiJ4iCl.C2.C3iC4»C5.C6i0ltl)2i03i04,05 
C3"M0N/0NE/AGtALtBB  ,SEtOE«  GK.     TTRP  ,PTRP  .DTRP .     TCRT  >  PCfiT,  DCRT 
CC*M0N/T(-RE/DPDT.D2P3T2.DLPDI_R.DTSDR.DTHDR,DPS1DR,      X6  1  ,  XB2  t  XO  1  .  XO  2 
Cai»MON/FI  \«E/     F  .T.OEN.DPDD.     E.H.S,     CV.CF.W.     WK 
COMMON/EIGhT/    DOLDT 
CONSTANTS     FFCM     MARShCAS     9/3/71, 
*RI T=<6, 99 ) 
}     FORMAT (30h        ***    CH4     ***       JANUARY     20     1S76  ) 

TTRP=90.68 
PT*J=0.  1  17«2Et75 
CTCP=2'3. 1472 
PCRT     =     45.956467  e*R 

TCRT=  190  .£; 
FCfiT=FS*TF (TCRT ) 
DCRT= 1 C. 15 
AC-=  .5 
Al_  =  .5 
&6=0.  8 
6E  =  4 .0 
Dc=l .2 

GK     =     0.083  1434 

CJ     =     100 

»K     =     101325/16.043 
A l=-4.  154Eje47 
A2=-4 .77  3901 66 
A3=3. 51 9S9C44 
A4=4. 257374  12 
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3  1=1. 74744€56 

32=2. 5767S£«5 

63=0.46538580 

Cl-S. 10335521 

C2=-li.420JJ74e 

C3=6.361 24273 

01=1.  £6363515 

02  =  -J. 54712477 

C3=l .a7406615 

RETURN 

£NC 

FUNCTION    PS*TF(T) 
METH4NE     VAPOR     PRESSURE     VIA     PRYDZ/GCCDWI N    DATA.     NOV.,     1570. 
NJTE.         PRESSURE     IN    BARS,  1.01225    BAR/ATM. 

IMPLICIT  REAL*8   (A-H.O-Z) 

C3MiaN/THRE/DPDT,D2PDT2.0LPDL.R.DTSDR,DTI-OR,DPSIDR,      XB1  ,  XB2.  XC1  ,X0  2 

GENERIC 

JATA        PTRP     •     0.1174356  75/.     TTRP/9C.68/.     TCRT/ 1 90.5 3/,      E/l.S/ 

OATA        A/4.7732553/,     3/1.7665875/,     C/- 0. £702812 /.    0/1.3311873/ 
I      XK=l-TTRP/TCRT 

*=(  1-77R°/T)/XK 

C=l-X 

IF<a>      2,3,4 
i     ^SATF     =      0 

CPCT     =    0 
RETLRN 

3  «  =  0 
Ml  =0 

GO     TO     5 

4  a    =    3**E 

II      =     -F*»/C 

5  DXJ)  r=7TRP/XK/T**2 

2  =  X*» 

Zl     =     *     +     >*" 1 
o     FZ     =     A*X     +     B*X**2     +    C* X**3     +    D*Z 
7     Fl     =      A    ♦     2*e*X     ♦    3*C*X»*2     ♦    D*Z1 
d     3iATF=PTRP»cXP(iTZ) 
JPCT=PSATF  4FHDXDT 
RETURN 
;NJ 
FUNCTION    OENLIQ(T) 
METHANE     SATD. LIQUID     DENSITIES     VIA     GCCDWIN. 
THIS    FUNCTION     REFITTED     VIA         .DENSATLQ.        4/29/71. 
IMPLICIT     FEAL*3     (A-t-.O-Z) 
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GENER  IC 

CGCMCN/EIGhT/     DCLOT 

l)ATA     TCRT/1S0.53/.DCRT/10. 15/ 

CAT*        A/O. £29403/.     T/l. 656635/.    C /-0.0 1 £387/.     E/0. 98/ 

1  X     =     T/TCBT 

2     =     1     -      > 
IF(Z .LT  .  .01  )     GO    TO     7 

2  F     =    -E*X**2/Z 

X3     =     EXP     (c) 

3  3FCX     =    -E*<2     ♦    X/ZMX/Z 

C     =     Z*»0.36 

4  G    =     A*Z     *     E*Q     +    O  X° 
OENLIC    =     CCRT* (1 +G) 

5  Gl     =     C*XP*DFPX    -    A     -     0.36*B*Q/Z 

6  CCLDT     =     DCPT*G1/TCRT 

RETURN 

7  IF(Z .LT.  1.E-0S)     Z=1.E-0S 
C=i*«  .36 

0=A*Z     +     9*C 

DENL 13    =     DCRT*( l.+G) 

Gl     =     -A     -     .36*E*Q/Z 

DJLDT     =     DCP7*G1/TCRT 

RETURN 
END 

FLNCTICN    TSATF(DEN) 
TSATOENI     FCfiMULA.         METHANE.     CONSTRAINED    TO    TRIPLE     PCINT.      IPTS- 
YIELJS     ALSC    ThE     FIRST    DERIVATIVE     RSP.    TL    RMODEN/DTRP. 
EON.,      (  7CRT/T-1  }  /AZ      =    F(R)      =     U(S)*<1     ♦     A1*LCG(R)     +      <  R- 1  )  *V»  (  R  )  )  , 
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C  ■*£«£•        U(S)     =     ( <S-l >/(ST-l > >**<8/3).        AZ    =     ( TCRT/TTRP- 1  >  t         AND     - 

Z  w(R)     =     A2    +     A2*R     +     ...     ♦     A6»R6. 

IMPLICIT     REAL*8     (A-H.O-Z) 

GENERIC 

COMMON/ONE/AG. AL.BB.BE.DE.  GK,     TTRP  ,PTRP  ,DTRP •     TCRT , PCRT, DCPT 

CCWM0N/THRE/0P0T,02PDT2,DLPDLR,DTS0R,DTHOR,DPSIDR,      XB  1  ,  XS2  •  XO  1  .  XD2 

01 NENSION    A (E) 

CATA     A     /    -0. £142449,     2.47S452S,     -e.759£>384.     41.05600105, 
1     -133. 4449220,     223.4955973,     -13  1.8267154,      58.e359o84/ 

1  AZ     =     TCRT/11RP    -     1 

CSOR    =     DTRP/DCRT 
SK     =    DSDR    -     1 

2  3=0EN/DTRP 
i=CEN/CCRT 

IF(  (  5-1  ) /SK.LT.O.  )     Q=-AES<  (S-t  i/SK  )**. 333333 
IF((S-I)/SK.GE.O.)     0=((S-  1  J /SK)**. 2 233  33 
IF IQ)     3,9  ,3 

3  u  =   o**a 

Ul     =     8*DSDR*0**5/SK/3 

4  V     =     ,7-1 
4=C 

*  1=0 

CO     6     K=2,3 

5  m    =     W     +     A(  K  )*R**(K-2> 

Ml      -     HI      -I     (K-2)*A(KI*R*»IK-3) 
a     CONTINUE 
3=0LOG(R ) 
F     -      U*(  1      4     A  (  1  )  *G     «-    V*W) 
7     Fl      =     01      +     A(1)*(U/R     +     U1*G)     +     U*V*H1     ♦     U*«(     ♦    U1*V*W 
3    TSATF=TCRT/( 1+AZ*F ) 

CT3  0P=-<A2/TCRT )*F1*TSATF**2 
kETJSN 
i     TSATF     =     TCST 
CT50F    =     0 

RETLRN 
END 
FUNCTICN    ThETAF(CEN) 
C  IF     S  1,     L     =    AG*<  S-l  )**3.  IF    5    1     It     U    s    -AL*  (S-l  )**3. 

C  YIELDS     ALSO     Tl-E    FIRST    DERIVATIVE     REP.     TO    RHODEN/DTRP. 

IMPLICIT     REAL*8      (A-h.O-Z) 
SEN5R  IC 

COMMON/ONE/ AG. AL.BB. BE, DE, GK,     TTRP , PTRP .DTRP •      TCRT .PCRT #DCPT 
C2M*ON/THRE/DPDT,D2PDT2.DLPDLR,DTSCR.OThDR.DPS  IDR,      XB 1 . XB 2 . XD  1, XD 2 
1     S=0cN/DCRT 
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DSJR     =     DTRF/DCRT 
IF( S )     9*9.2 
!     C=S-l 

J2=Q**2 
G3  = J**3 
1F(0)     3,8.4 
)     u     =    AG*Q3 

Ul     =     3»AG*C2*CSCR 
SJ    TC     5 
t    J    =     -AL*Q3 

Ul     =     -3*i»L*C2*DSDR 
j     IF10.L7.-5C.)      L=-50. 
XP    =     EX°(U  ) 
TS     =     TSATMCEN) 
THETAF     =     1J*XP 
'     CTI03     =     TS*U1*XP    +     0TSDR*XP 

RETURN 

I     ThETAF     =    TCRT 

CTHDK    =     0 

RETURN 

I     ThETAF    -     0 

07H3R=l.E     05 

RETURN 

Ero 

FUNCTION     PHIF(T) 
X3     =     PHI        =        X*(  l-EXP(-U)  )»         U    =    BE     +    BE/X, 
YIELDS     ALSC        CPhl/OR.     DPHI/OX,     AND     02PHI/DX2. 

IMPLICIT    REAL*e     (A-M.O-Z) 
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GEr>  =  RlC 

CG*MON/CNE/AG.AL.SB  .BE.OE.GK,     TTfiP  .  PTRP .DTRP.     TCRT , PCRT . CCRT 
CJMM0N/THPE/0PDT.D2°DT2,DLPDLR.D1 SDfi.D THDR.DPSIDR.     XB1.XB2.XC1.XD2 
X=T/TCRT 
U=BS+BE/X 
U1=-BE/X<*2 
o2=-2*Ul /X 
I     XP=EXP     (-U) 
2=1-XP 
Zl =U1*XP 
Zi=<U2-U  1**2)*X? 
1     PHIF     =     X*2 

X81      =      X*Z  1     +     Z 
XB2     =     X*Z2     +     2*Z1 
I     RETURN 
END 
FUNCTION     PS  IF (T. DEN) 
XJ>     =    PSI         =         (1-W*LQG<1  +1/W)  )/X,  h     =     (T-TH)/'(DE*TCRT  ). 

YIELDS     ALSC     XO  1     >    DPSI/DX,      XD2     >    D2PSI/DX2.     AND    OPSI/OB. 
IMPLICIT     RE*L*8     (A-t-.Q-Z) 
GENER  IC 

CCMMON/ONE/AG. AL.B3.BE.DE. GK,     TTRP .PTRP  .DTRP t     TCRT , PCRT . D CPT 
CC*vlON/ThRE/C°DT,D2PDT2,OLPDLR.DTSDfi.DTHDR.DPSIDR.      X3  I  .  XB  2  .  XD  I  •  XD2 
1      TH     =      THETAF(DEN) 

A     =     (T-Th J/DE/TCRT 
IF{W)     2. 2,  3 
I     °SIF=1 
XC1=0 
XO2=0 
DPSIDR=0 
PETUPN 
J     x=7/TCPT 

D#CR=-OTHDP/DE/TCRT 
D*CX  = 1 /CE 
»     U=l/X 

C'JCX=-U/  X 
l)2JDX2  =  -2*CUCX/X 
G=DLOG(  1  .  ODC+1  .0D0/\») 
«3     =      1+W 
V      =      1-V»«G 

»   jvd*    =    l/ws-i 

C2VO*2     =      l/'*/WS**2 

PSIF     =     L*\ 
'     DPSIDR     =     U*DVDW«DWDR 
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xO I     =     U»D\D»*D*DX     +     V*DUDX 
3     XC2    =     U*02VC*2*D*DX»*2    +     2*DUDX*Q  VD«)*U  *DX    +     V*J2UDX2 
9     KETURN 
END 
FUNCTION     DFDRF(T.DEN) 
NOTE.        DPDRF     =    PRES5UPE.     EAR. 

REDUCED     DERIVATIVE.        DLPDLR     >    R*(DP/DR)/P        IS     IN    COMMCN. 
(2-t»*X/R    =     F(T,RJ     YIELOS     ( DP/DR » /OTRP/GK/T    =     l+(R2*Fl+2*S*F )/X. 
IMPLICIT    REAL*8     (A-H.O-Z) 
GENE^  IC 

CCKMON     Alt  A2.  «3i A4 . 31. B2. B3. B4 1 C 1 • C2. C3.C4, CSt C6.D1.D2.D3.D4.D5 
COM^DN/ONE/AG.AL.BB  .9E .DE.GK,     TTRF  ,  PT.JP  tDTRPt     TCRT  ,  PCKT  ,  DCRT 
CJMMON/THRE/DPDT.J2POT  2. DLPOLR . DT SDR. DTHDR. OPS IDR .      XB 1  , XB 2  .XD 1  .  XD2 
CQMrtuN/FCUF/FPT.DFDX.DFDRl 

1  X=T/TCRT 
S=CEN/DCRT 

DSDR     =     OTPF/DCBT 

2  R=CEN/DTRP 
fi2=R*»2 

3  XE    =     PMIFIT ) 

XC     =    1  /x 

4  XD    s     P£  Ic(  T.DEN) 

XDO    =     CPS  1CR 

6  A     =    Al     +     A2*R    +     A3*R2     ♦     A4*R3 

AC     -     A2     +     J«A3*R     +     3*A4*R2 

7  e     =     81     +     E2*Fi     +     E3*R2 

33     =    32     +     2*B3*R 
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3     C    =     Cl*fi     ♦     C2*R2     ♦     C3*R3 

CD    =    C  I     ♦     2*C2*R    +    3*C3*R2 
1  a    U    =    S-l 

LI      =    OSOR 
il     V    =    01     +    D£*R     ♦    D3*R2 

VI     s    Q2     4    2*C3*R 
14    3     =     U*  V 

JC     =     U*V1      ♦     U1*V 
F     =     A     +     E*XE     ♦    C*XC     +     D*XD 

16  OFJR     =     AT     +    SD*XB     +     CO  *XC        +       D*XDC     ♦     OC*XO 

17  G    =     1      ♦     (R2*DFDR     ♦     2*R*F)/X 

G=l +F»R/X 
OL3QLR     =     C/O 
Z 0     SFCRF    =     CEN»CK*T*Q 
DFOR1 =  DFOR 
RETURN 
END 

END     JF    FIl_E     CN     S>£INi 
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C  rIEGRESSICN     FRCGRAM     TO     FIND    T*.V*.      AND    K(IiJ)     USING    PURE 

C  COM33NENT    DATA 

IMPLICIT     REAL*8     (A-M.O-Z) 

OIMSNSION    DIM  (3)  .DIN2  (8).DIN3(8).NC(20).T(20,  110),P(20.110>. 

*  ^0(20.1  lOt.RKTI  20.  110).ANU(20t5)  ,3NU( 5). TS(S)  .  VS (5 )  . AK (5  ,5  )  • 

*  NCJR<6.5).F(2000), X{ I  5  ). W { 1200 Oi . CT( 20  >. CV( 20 ) .RE S 1 ( 2  0. I  1  C ) . 
«    RE32I20.110) ,NA,ME<25) . NTI TLE ( 20 ,25 ) 

C  LIMITATIONS 20    COMPONENTS 

C  5     GRCJP3 

C  110    DATA     POINTS    PER     CCKFONENT 

C  UNITS PRESSURE     IN     ATMS    OR     EARS     (OPTION     IB) 

C  VOLUMES     IN     CC/GMCLE 

C  TEMPERATURES     IN    KELVIN 

C  NCTE3 NAME     IS    FORTY     CHARACTERS 

C  LSE    ONE    BLANK     CARD    AT    THE     END     OF     A    COMPONENT     CATA     SET 

C  USE    ANOTHER     BLANK    CARD    AT     THE     END    OF     ALL     THE     CATA 

C 

C  SET     JP    COMMCN    BLOCKS 

C3MM0N     /DA7A/     T. P. RO .RKT. ANU 

CCMMCN    /PAPAM/    TS.VS.AK 

COMMON     /DECIPH/    NCOR.NG »ND ,NNM 

COMMON    /RESULT/    RES1.RES2 
C  NOPT    DETERMINES     IF     A    REGRESSION    OR     A     CALCULATION     IS    TO     BE    PERFORMED 

C  IF     JOPT     IS     0     THEN     THE     REGRESSION     PROCEEDS    ANY     CTHE  R    VALUE 

C  CF     NOPT     CAUSES     ONLY     A     CALCULATION     TO    3E    PERFORMED 

*EAD( 5.47) KCPT 
47  F03MAT(  I   1  ) 

C  FIRST     INPUT     IS     THE     NUMBER     OF    DIFFERENT     GROUPS 

Sitf    READf 5. 1 .END=9999> NG 

1  FO  fi  M  A  T  (  I  2  ) 

Z  INITIALIZE     COUNTER     FCR     NUMBER     CF    MOLECULES 

N*M     =      0 
C  ENTER     DATA     FCF     FIRST     COMPONENT 

2000        REAJH5.2)      hAME.TC.VC 

2  FCfi'-IAT  C25A*.  2F  10.0  ) 

C  CHECK     TO     SEE      IF     ALL     CATA     HAS     BEEN     INPUT 

IF(TC  .LT.l  .ODO)     GO     TO     7777 
C  NJ     IS     THE     STC ICHIOMETR IC    VECTOR    FOR     THE    COMPONENT 

READ! 5.2  5) (ENU U ) . I =1 . NG) 
25  FORMAT(  5F  ICO) 

C  UPDATE     CCKFCNENT     COUNTER 

NNM     =     NNM    +     t 
C  SET     COUNTER    FCR    NUMBER     OF     DATA     POINTS    FCR     EACH     MOLECULE 

ND(NN^)      =     0 
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C  JEGIN     fVCKlNG     INPUT    3A TA 

CC    50     1=1.10 
53  .NTITLECNNM  .  I  )     =    NAME(I) 

CT(NNM)     =     TC 

Cv(NNK)      =     VC 

JD    75     1=1. NG 

75  avj(nnm.i)    =    enu(I) 

z  cntsr  optio  fcr  pressure  units 

REA0(5.3)      IE 
3  r  C  S  M  A  T  (  1 1  ) 

100J        REATX 5. 1)      TEMP 
i  FCftMAT(F10  .0 ) 

C  CHECK     TO    SEE     IF    END    CF     A     CCMPCNENT     CATA    SET 

IF(TEMP.LT.  1.0DO)     GO     TO     2000 

REAJ(5,5)CIM.DIN2.DIN3 
5  FGHMATf 8F 1 C. 0) 

C  ChECK     TO     SEE     IF    TOO     CLOSE     TO     CRITICAL    TEMPERATURE 

ChECK     =     DAEStTEMP/TC    -     1.000) 

IF(CHECK.LT.O.OSDO)     GO     TO     1000 
C  ChECK     ON     PRESSURE    UNITS    AN!)    CONVERT    3ARS     TO    ATMS 

IF(IB.EQ.O)     GC    TC    3000 

jo    too    1  =  1. e 

100  OIM(I)     =     CIMCD/l  .01325 

C  CHECK    FOR     2EPO    DATA    POINTS.     CCUNT     AND     PACK 

3000         CC     200     1=1.8 

IFIOIM  (  I  )  .LT.l  .000  .OR  .DIN2(  I  )  .LT  •  l.OUO  )     GO     TO     200 

IFC VC/DrN2( I ).LT.l .300)     GO     TC    200 

NO(NNM)     =     NO(NNM)      +     1 
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NOL*     =     ND(NNM) 

T(NNM,NDU«)     =     TEMP 

PINNM.NDUK1      =     DIM   (I) 

RO(NNM.NDUC)     =      I . OD 0/0  I N2 (  I  ) 
?30  CONTINUE 

Z  RETURN    FOR     ANCTHErt     ISOTHERM 

GO    TU      1000 
C  AFTER     ALL     CATA     IS     IN    ENTER     KNCWN    PARAMETERS     AND    FLAGS     FOR 

C  THOSE     TO    BE    FOUND     BY    THE    ROUTINE 

C  ENTER     A    NEGATIVE    VALJE     OF     T*    OR     V*     IF     IT    IS     TO    BE     FOUND     BY     THE 

C  PROGRAM.         THE     INITIAL     GUESS     WILL     EE    NUGAT IVE    THE     INPUT     VALUE 

C  ENTER     K>1     IF     IT     IS     TO     BE    FOUND 

C  FINO     TOTAL     NUMBER    OF    DATA     POINTS 

7  777        NSUM     =     0 

CO     300     1=1. NNM 
300  NSLM     =     NSU*     ♦     NO ( I ) 

C  ENIER     T*    AND     V* 

RcAJ(5.6)(TS(  I).I=1.NG) 

REAO( 5.6 ) ( VS( I) .1=1 »NG) 
o  FOR>4AT(  5F10.0) 

C  FINO     NUMeEF     OF    UNKNOWNS     AND    SET    UP    COf-.RELATI  ON    MATRIX    TO     BE 

Z  ABLE     TO     RETRIEVE    THEM     LATER 

C  INITIALIZE     COUNTER     FOR     NUMBER     OF    REGRESSION    PARAMETERS 

1987        NF  =  0 

NGP     =     NG     +     1 

CC    400      I=1,NGP 

GO    4  0  0     J=l .NG 
400  NCORC I.J )     =     0 

DC    500      [=1.NG 

IF(TS(I ) .GT.l . ODO)      GC     TO    5000 

NP     :      N°      ♦       1 

NCCR( 1 . I )      =     NF 
C  SET     INITIAL    GLESSES    AS     *Z     GC 

X(N=>  )     =     -I  .0O0*TS(  I  ) 
500J        IF( Vj< I ) .GT.l .ODO )     GC    TO    500 

NP      =      NP      t       1 

NCC^( 2.1)     =     NF 

X( NP)      =     -1 .0CO*VS< I ) 
500  CONTINUE 

C  NCW     WILL     F.EAC    K'S     AND     FINO    UNKNOWN    ONES    AND     ENTES      INITIAL     GUESSE 

C  NOTICE     THAT    ONLY     THE    FIRST     NG- I      RCkS     OF     AK     MUST     BE     ENTEREC 

C  REMEMBER     TO     ENTER     K> 1     IF    UNKNOWN 

AK( 1.1)      =     C.ODO 

IF(NG.EQ.l)     GC     TO    3333 
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NGX      =     NG     -     1 

30     60  0     I=1,NGM 
>00  KEA0I5.7  )  (AK(  I.  J).  J=U  NG) 

7  rCKMAT(5FlC.C) 

CO    650     1=1. NG 
650  AK(NG.I)     =     0.0D0 

JO     70C     1=1 iNGM 

IP*I+1 

DC     700     J=IP,NG 

IF(  AK  (  I.J). LT.  1.000)     GO    TO     700 

NF     =     KP+1 

X{NP)      =     O.CDO 

IP2    =      1+2 

NCCRCIP2.J)     a     NP 
7  00  CONTINUE 

C  CI-ECK    TO    SEE     IF     IN     REGRESSION     OR     CALCULATION    MODE 

333J         IF( NJPT. NE.O)     GC    TC    3335 
C  3±GIN     REGRESSION 

CALL     VA0  5AC(NSUM,NP,F.  X.  1.0D-1.I.0D2.  I .00-3.99.  10.  W ) 
C  *FT=R     REGRESSION    COMPLETE     CALCULATE    FUNCTIONS     ONE     MORE    TIME    FOR    GUTF 

J335        CALL     CALFUMNSUM.NP.F.X) 

C  CALCULATE     AVERACE    ERROR    FOR     ALL     DATA    PUINTS     REALIZING     THAT     THE 

C  FIRST    DATA     POINT    ON    EACH     ISOTHERM    HAS    ZERO     ERROR     ANC    SHOULD    NOT     8E 

C  CCUNTED 

NZ     =     0 

SLIE^R    =     O.OOC 

CC    300     I=1.NSUM 

St-MERR     =    SCVERP    +    DA8S(F(I)) 
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IF(F(  I  )  .NE.O.ODO )     GO    TO     300 

NZ     =     NZ+-1 
330  CONTINUE 

AVGABS     =    SUMERR/CNSUM-NZ) 
C  3EGIN    GENEP-AL     CLTPU7 

«H ITE( o,    1  1) 

11  FCPMAT (1  HI  .20X, • SUMMARY     OF     RESULTS     •///  ) 
MR IT5( 6. 12  > 

12  FCfi«tAT  (20X,  "CHARACTERISTIC    PARAMETERS*/// 

»     20X. 'GROUP     NUMBER"  .10X,  •     TSTAR     ",10X»"     VSTAR      •// 
MRITE(6.  13  )(  I.T5<  I  ),VS(I).I=i.NG) 

13  FCPMAT  (2  5X,  12,  16X, F7.2 . 1  OX. F7 .2 ) 
*RITE(6. 14) 

14  F3fiMAT<///. 20X, •     K        MATRIX     •///) 
OC     1200     1=1 ,NG 

1200        XRITSt 6.  15)(AK{ I, J)  .J=l  ,NG) 

15  F3RMAT(20X,5(F10.4. 5X) .///) 
MNN     =     NSUM    -     NZ 

MR  ITE(6t   16) AVGABS. NNN 

16  FOPMAT (20X , "THE     AVERAGE    ABSOLUTE     ERROR     MAS     ",F15. 

*  20X,      "CALCULATED    FOR     '  .13.     •      CATA    POINTS     • /// ) 
C                 3£GIN     INOIVIDUAL    COMPONENT    PRINTOUT 

NF=0 

DO     »0C     11=1. NNM 

•RITE(6,  17  J 
1  7  FGRMATUHl  ) 

MR  ITi(6,  18)(NTITLE{  II .  I  ) . I = 1 .25 ) 
13  FCPMAT < 10X ,25A2,/// ) 

»RITi(6.19)CT(I£l.CV(II) 
H  FOR"*AT(  10X,  "CRITICAL    CONSTANTS"// 

*  10X, • TEMPERATURE •  ,  F7 .2/ 

*  10X, • VOLLM" «,F7.2//) 

MR  IT  E(  6,  2  1  )(ANU(  II,  I).  1  =  1.  NG) 

21  FCSMATC I  OX  ."STOICHI  CMETRI C  VECTOR".//. 

*  15X, 5F7.4///) 
MKITE (6, 22 ) 

22  FCRMATOOX.'TEMP'^X.1  PRESSURE"  .2  X  .  •  VOLUME"  .  2X  ,  •  C 

*  "0EL°X",4X,  'DEL  PC  • . 4X.  "ERROR" . 4X • • SERKOR •  ,///) 
NOLM=ND( II ) 

SbMI=  C.  00  0 

NCCJNT     =     0 

DC      1300     LBjM  .NUOM 

>4F  =  NF  +  1 

VCL     =     1 .OCO/RO ( II.LEJ ) 
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NNC-M  =  NF 
13500     cfn     =     R5S2  (  IItLBJ)-RES 1  (  I  I  tLBJ  ) 

aR{TE(6.23)T(  II.L3 J)  tP( II  tLBJ)  .VCL.fiO( I I.LBJ>t RES!  (II.LBJ  ). 

*  RE32(  I  IiLEJ  ),ERR.F(NND'JM) 

2  3  F0RMAT(9X,F7.2.2X,F7.2.1X.F7.3.2X.F7.4,lX.F8.3.1X,F8.3.1X.Fa.3. 

*  IXiFS.3) 

SUMI  =  SUM  I  +  C»ES (FINNDUM  ) ) 

IF(F( NNOUM) .NE. 0.000)      NCOUNT     =     NCCUNT  +  1 
1300         CONTINUE 

AES^I     =     SUVI/NCCUNT 

«RIT5( 6. 24)AERRI 
24  FUK4AT(///, 20X, 'AVERAGE    ABSOLUTE    ERROR     =     • . F 1 0. 4 ,/// ) 

900  CuNTI NUE 

GO    71     995 
9993     STOP 

£ND 

SUBROUTINE     CALFJNC NSUM.NP.F.X) 
C  THIS     SUBROUTINE     IS     USED     BY     VA05AD     FOR     THE     REGRESSION 

C  ALL     UNITS    AND     RESTRICTIONS     ARE    THE    SAME    AS     FOR    THE     MAIN     ROUTINE 

IMPLICIT     REAL*e     (A-t-.O-Z) 

JIV£NSION     NC(20).T(20.110).P(20,110).RC(20.110).RKT(2G.  1  I C ). 

*  ANJ(  20. 5)  ,  TS( 5)  .VS(5)   .AK15.5)  . NCCR(6.5).F(2000).X(15). 

*  3E.S1  (20.  1  10  ).RES2(20.   110>.ANUI(5).TI(110).PI(  110)  ,  RO  1(110) 
C                  SET     UP    COMKCN    BLOCKS    FOR    COMMUNICATION 

COMMON    /DATA/     T, P. RC .R K T . A NL 
CCKXCN    /PAP/M/TS.VS.AK 
COK.HCN    /DECIPH/     tsCCK.NG.ND  ,NNM 
C3MON     /RESULT/    RES1.RES2 
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CC*MON     /SAVE/     III 
C  FIRST     SORT     X    VECTOR     TO    RECOVER     PARAMETERS 

1 963        JC     100     J=l .NG 

IF(NCOR(  1. J) .NE. 0)      TS(  J)=X<NCCP(1  .J)) 

IF(NCOR{ 2, J  >.NE.O)     VS(J)=X(NC0R(2.J)) 

1F( NG.EQ.l )     GC    TC     1 00 

NGP     =     N3     +      1 

CC     110      1=3, NGP 

I *  =  I-2 

IF(  4COR(  1 1  J).NE.O)      AK(  IM, J)=X(NCOR<  I. J)  J/10O.OO0 
1  10  CONTINUE 

100  CCNTI NUE 

IF(  MG.EQ.l  )     GO    TO     250 

NG»>     =     NG-1 

DC     200     I =1  ,NGM 

IP    ~      I  +  l 

JO    200     J=IP,NG 
200  AK< J,  I  )     =     AK(  I  .J) 

C  CONTINUE    TC    RECOVER     PARAMETERS     WHILE    CALCULATING    FUNCTIONS 

C  GLTER    LCCF     PEF.FGRMS     CALCULATIONS     FOR    ALL    MOLECULES 

250  HF     =     0 

OC     300     11=1. NNM 

NOLM=ND(  I  I  ) 
Z  SuRT     STOICHIOMETRIC     MATRIX 

JC    tOO     JJ=1.NG 
400  ANLI  (  J  J)  =AM(  I  I  •  JJ  ) 

C  SORT     T.",RhO 

JC     500     KK= 1  .NCUM 

TKKK  >     =     T(  I  I  ,KK) 

P  IKK  }     =     P(  I  I.KK  ) 
500  RGI(KK)      =     PC(II.KK) 

C  TOLJ     AND     III     ARE    USED     TO     PREVENT     UNNECESSARY     CALCULATIONS     FCR 

C  TcfiHS     THAT     AUt    ONLY     TEMPERATURE     DEPENDENT     IN    THE     EOS 

TCL3  =  C. ODC 

DC     JOO    L=  l.NDUM 

|\F  =  NF  +  t 

IFOAHStTI  (L)-TOLO  )  .L.T.0.01D0)      GO     TO    6000 

I  I  1  =  0 

TCLO=TI (L) 

TR  =  T I  <  L  ) 

PR  =  o I (L  ) 

RfisRO  I  (L  ) 

CALL     GP(  NG .RR.ANUI  ,  AK.  TS , V S. TR , PRC i 
6300  111=1 
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TM=TI (L) 

PM  =  ^>I  (L  ) 

fi*=fiul (L) 

C4LL     GP(NG.CM,ANUl.AK, T S . V S .T K . PMC ) 

RES  1(  ILL  )     =     PM-PR 

RES2 (Ii.LJ     =     PMC-PRC 

NMD  UM=NF 
6750        IFCRESi  (  ILL  )  .GT  .1  .000  )     GO     TO     C760 

F(NNi)UM)      =     O.ODO 

GO  TO  200 
6760   F(NNDUM)  =  (  P ES 2 ( I  I • L) /RES !( I  I .  L  )  -  1 . ODO ) *  100. ODO 

200         CONTINUE 

RETURN 

ENC 

SLBROLTINE  GP <NG .RHC. ANU . AK ,TS • VS «T .P ) 

c     pjro-ise:   to  calculate  the  pressure  of  a  pure  component 

C  LSING    THE     GROUP    CCNTR  IBUT ION     EQUATION    OF     STATE 

c  variables:       NG-THE     NUMBER    OF    DIFFERENT     types    OF    GROUPS 

C  RHO-THE    density    in    GMOLE/CC 

C  ANU-THE    VECTOR     CF    STOICHIOMETRIC     COEFFICIENTS 

C  AK-THE     MATRIX    OF    BINARY     INTERACTIONS    CONSTANTS 

C  TS-THE    VECTOR    OF     CHARACTERISTIC    TEMPERATURES     IN    K 

C  VS-THE      VECTOR     OF     CHARACTERISTIC     VOLUMES     IN     CC/GMCLE 

C  T-TKE    TEMPERATURE     IN    K 

C  P-THE     PRESSURE     CALCULATED     IN     ATM 

C  SUBROUTINES    REQUIRED:         SIGMA 

C  VIRIAL 

C  BPHS 


145 


UGUST     26.      I9ei  11126  NEKDC       SYSTEM    SUPPORT     UTILITIES 


IMPLICIT    REAL*9     (A-H.C-Z) 

DIMENSION     R0<  10),ANU(5),AK(5.5).TS(5)  iVS(  5)  , B(  5,5)  ,SIG(5)  ■ 
*     SIoMllO, 10 ) 
C  RG     IS     THE     \ECTO=    OF     GRCUP    DENSITIES 

C  3     13     THE    MATRIX    OF     GENERALIZED     VIRIAL     COEFFICIENTS 

C  SI3     IS     THE     VECTCR     OF    HARD     SPHERE     LIAMETERS 

C  SIGM      IS     THE    MATRIX     OF     HARD     SPHERE     DIAMETERS 

CCVMCN     /SAVE/      III 
C  III     IS     USED     TC    DETERMINE     IF    THE    TEKPErtATURE     HAS     CHANGED 

C  IF     III      IS     C    ALL     CALCULATIONS    ARE    PERFORMED 

C  IF     III     15     NCNZERO    ONLY     Tl-E     HARD    SPHERE    PRESSURE     IS     CALCULATED 

C  HI     SHOULD    9F     0    FOR     THE    FIRST    CALL    CN    AN     ISOTHERM 

C  CCNSTANTS     ARE     DEFINED 

PI =  3. 141 5926525D0 

R=82.0S6D0 
C  CALCJLATE    THE    TCTAL     NUMBER     OF     GROUPS    AND     THE    GROUP     DENSITIES 

TNG  =  O.ODO 

DC     10     1=1. NO 

TNG     =     TNG     ♦     ANU(I) 
10  SGI  I  )      =     4NU  I  )*RHO 

C  CHEC<     TO     SEE     IF    TEMPERATURE    HAS    CHANGED 

IF{ II  I .GT.O)     GC    TO     100 
C  CALCULATE     HARD     SPHERE     DIAMETERS    AND     VIRIAL    COEFFICIENTS 

CALL     SIGMAING.T.TS.  VS.  SIG) 

CALL   VI RI AL (NG.T .TS. VS. AK, B> 
C       CALCJLATE  PEPTUP.BAT  ION  TERM 

VC  =  0. ODO 

DO  30  1=1. KG 

CO     30     J=1.NG 
30  VC     =      VC     +     ANUII )* ANU( J )*B( I • J ) 

C  CALCULATE     THE     HA  RD      SPHERE     CONTRIB  L  T  I0.-( 

100  CALL     eFHS  INC. RG.S I G. HSBP ) 

C  ADD     ALL     CCKTRIBLTI CNS     TOGETHER 

BP  =  HSeP-PHC*PHO*VC/2.0D0  +  PHOi'(  1  .ODC-TNG  ) 
C  CALCJLATE     THE     ACTUAL    PRESSURE 

3=a°*R*  t 

C  SET     III>0     FGR     NEXT     °*SS    THROUGH    THE    ROUTINE 

I  I  1=1 

RETURN 

ENC 

SLEPCUTINE    SIGMA( NG.T.TS.VS.SIG) 
C  THIS    SUBROUTINE     WILL    CALCULATE     THE     VECTCR     OF    HARD     SPHERE 

C  DIAMETERS     USEC     IN    ThE     HARD    SPHERE    PRESSURE    CALCULATION 
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C  THE    FUNCTIONAL     FORM     USED    FDR     THE     TEMPERATURE    DEPENDENCE 

C  15     A     WCDIFICATIGN    OF    THAT    °RGPGSED     BY     BIENKOWSKI     AND 

C  CHAD     FOR     THE     INERT     GASES 

C  VARIABLE    L  1ST: 

C  NG-ThE    NUMBER    OF    GROUPS 

C  T    -THE     TEMPERATURE     IN    KELVIN 

C  TS-THE    VECTOR    OF     CHARACTERISTIC     TEMPERATURES 

C  FOR    THE     GRGUPS     CONSIDERED     IN    KELVIN 

C  VS-THE     VECTOR    OF     CHARACTERISTIC     VOLUMES     FCR 

C  THE    GROUPS     IN     CC/GMOLE 

C  SIG-THE     VECTOR     CF     HARD    SPHERE     DIAMETERS     THE 

C  UNITS    ARE     SUCH     THAT    KHG*SIG**3     IS 

C  CIMENSICNLESS     IF     RHO      IS     IN     GMOLE/CC 

C  THESE     ARE     NO    CTH5R     SUBROUTINES     REGUIKED 

C  LIMITED    TO     5     GROUPS    THIS    CAN     EE    CHANGED    BY    ALTERING 

C  THE     DIMENSION     STATEMENT 

IMPLICIT     REAL*8     (A-H.O-Z) 
C  ALL     CALCULATIONS     IN     DOUBLE    PRECISION 

DICiENSION     1S<5)  i  VS(5)  .SIG(5) 

DIMENSION     XC(6) 

CAT  A     XC/O.  I  36  32 DO. 0.2 644  6 D- It  0.80  0  57D- 1.-0.899  99D- I. 0.390  6SD-1 
*       -0.609840-2/ 
C  DEFINE     CONSTANTS     REQUIRED 

PI     =     3. 1415S26535D0 

T3    =     1.0D0/2.0D0 
C  ALL     CALCULATIONS     IN     ONE    LOOP 

OQ     10     I =1 .NG 

TR     =      T/TSI   I  ) 
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20  TRI=l.CDO/Tfi 

F=XC(1)  +  TR1*(XC<2)+TRI*(XC(3>  +  TRI*(XC<4)«-TRI*(XC{5)+TRI*XC<6>)))) 
10  5IG(X  )= (6.0CO*F*VS(  D/DI  )**T3 

RETURN 

ENC 

SUBROUTINE     VIRIALCNG.T .TSi VS. AK , B ) 
C  PURPOSE:         TC     CALCULATE     THE     MATRIX    OF     VIRIAL     COEFFICIENTS 

C  VARIAeLES:         NG-NJMBER     OF     DIFFERENT     TYPES    OF    GROUPS 

C  T-THE     TEMPERATURE     IN     K 

C  TS-THE     VECTOR     OF    CHARACTERISTIC     TEMPERATURES 

C  VS-ThE     VECTOR    OF    CHARACTERISTIC     VOLUMES     IN    CC/GMOLE 

C  AK-THE     MATRIX     OF    BINARY    CONSTANTS 

C  B-ThE     MATRIX    OF    VIRIAL     COEFFICIENTS     IN     CC/GMOLE 

IMPLICIT     REAL*8     <A-H,0-2) 

JIMfNSION     TSI 5) . VSC5). AK(S.S)  ,B ( 5 .5) . TR( I  0 . 1 0)  . VR( 10.  10) . 
S     Alb) 
C  TR     I S     A     MATRIX    OF     CHARACTERISTIC     TEMPERATURES 

C  VR     IS     A     MATRIX    OF    CHARACTERISTIC     VOLUMES 

C  A     IS     t     VECTOR     OF    CONSTANTS     FOR    THE    VIRIAL     COEFFICIENT     CORRELATION 

J  AT  A     Ay 5. 0  65  ID 0,-2. 32 101. 5. C8270 1  .-*. 0662D1 . 1. 5fcl BQ1 .-2.3 1 7  300/ 
C  DEFINE     ONE     CONSTANT     NEEDED 

T3=l. 0D0/3.CDC 
C  ALL    CALCULATIONS     INSIDE     ONE    LOOP 

OC     10      1=1 . KG 

DO      10     J=I.NG 

IF(  I-J  >      1.2.1 

1  TMI,J)=OSCRT<TS(I  )*T5  (J  )  )  *  (1  .ODO-AKl  I,  J)  ) 
VR< I . J>=( ( ^£(I)**T3+VS( J)**T3)**3)/8.0D0 
GC    TO     3 

2  TR(I . J) =TS( I ) 
VR( I . J)=VS(  I  ) 

3  TI=TR  (I,  J)/T 

d(I.J)  =  (A(l)  +  TI*tA(2)+TI*(A(3)+TI*(A(4.)+TI*(A(5)+TI#A(6))))))* 
*        VR(  I.J  ) 

e ( j. i )  =  e<  i. j ) 

10  CONTl NUE 

RETURN 

ENC 

SUBROUTINE     BPH S ( NG , KG . SI G . HSBP ) 
C  PURPOSE:         TO     CALCULATE     THE    PRESSURE    OF     A    COLLECTION    OF     HARD 

C  SPHERES     USING    THE    GENERALIZED       HARD    SPHERE     EQUATION 

c  mote:     THE    CLTPLT    is    actually   beta*f    IN    GMCLE/CC 

c  variables:      ng-the    number   of  different    types  of  groups 

C  RG-VECTCR     OF     GROUP     DENSITIES     IN     GMOLE/CC 
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C  SIG-THE    VECTOR     OF     HARO     SPHERE     DIAMETERS 

C  HSBP-DETA*P    FOR    THE    HARD    SPhERES     IN    GMOLE/CC 

IMPLICIT     REAL*6     (A-H.O-Z) 

OIXEN3ION     PG<1  0  >tS  1C(5  ) 

DATA     AL/-4.2D0y' 
C  ENTER     THE    VALUE    OF    PI 

Pl  =  3.  1*1 592652500 
C  ZSR3     REGISTERS    FOR     CALCULATION     CF     SUMMATIONS 

20  =  0 .OCO 

zi-a.  coo 

Z2=0.0C0 

Z3= J. ODO 
C  CALCULATE     SUMS 

DC     10      I = 1 1  N  C 

AD0=PI/6. OCO*SG(I) 

ZC  =  ZO-fACD 

Z1=Z1 +ACD*S IC ( I ) 

Z2=Z2+A0D*SIG<  I)*3IG(I  ) 
1 0  Z3=Z3+A0D*S  IG(  I  >*S  IG(  I  )*SIG( I  ) 

3UM  =  l  .0Q0/(  1.0D0-Z3) 

P  1  =  Z0*DUM 

P2  =  J .OCO  *Z  1*Z2*CUM*0UM 

P^=Z2*Z2*Z2*DUM*DUM*CUM* (3 . 0  00+ AL* Z3 ) 

HS3P= (6.00C/PI )*(P 1+P2+P3) 

hETURN 

END 

ENU    JF     FILE     ON     SYSIN. 
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C  PROGRAM     TC     CALCULATE    REDUCING     VCL.UXE    CORRECTION    FOR    ThE    H#FD 

C  CONVEX    BODY    EQUATION    OF     STATE 

C  VEfiSICN     CF    2MABCH     1980 

C  INPUT     VARIABLES    ARE      : 

C  1  .     AL  SHA°E     FACTOR    OF     INTEREST     (fi*S/3V) 

C  2.     a        LOHEfi    LIKIT     OF     INTEGRATION     FCR    PACKING    FRACTIOIN 

C  3.     T     UPPER      INTEGRATION     LIMIT 

C 

C  SLEKOLTINES     AND    FUNCTICNS     PEQUIREC 

C  RacADC.HVALFUNCTN.FV.G V.GALSS 

IWFLICIT     REAL*8      (A-h,0-Z) 

DIMENSION     FCL(  1  0)  .  I  NA(3) 

COMMON     X, AL 

EXTERNAL     FLNCTN 
1000    REA .3{  5.  1  «END  =  S0O  )     AL 

1  FGR-IAT  (DIO  .3  ) 
3=C.2D0 
T=0.S5O0 

X=0 .9D0+0.  144C0*AL 
IT£ST=2 

10  CALL     FQUADC (  ITEST,  X,F,   100.  1.0D-S. l.CD-5.0. 2D0.R0L. I NA) 
GO     TO     (1  1  .  12  .13  ,14  )  .  ITEST 

11  CALL     5VAL(   > .F ,AL .B . T) 

ec   to    10 

12  »RITE<6,2)     X.F.AL.B.T 

2  FJ3RMA7<  1H0  .  10X.  'BE  ST    VALUE     CF     A=«,C15.5// 

*  10X.»     THE     FUNCTIONAL     AT    ThE    MINIMUM     IS    =«.D15.5// 

*  10X,"     THE     SHAPE    FACTCR     IS     ='.D15.5// 

*  10X.    'THE      INTEGRAL      »AS     EVALUATED    F RCM* .Dl 5. 5 . • TO « . D 15.5// /// ) 
iC    TC    1 000 

12     «RI  TEC  6.  2)  >.F 

3  F0HMAT(10X,'     ERROR     OUE     TO     ROUNDOFF        A=  •  ,0 1 5. 5.  »F=»  ,D 1 5. 5/// // ) 
GC     TC     1000 

14     ■SITS(6.4)X«F 

4  FCfiMAT (  10X ,  'TCC    MANY     ITERATIONS        A=  •  . D  1  5. 5.  »F= • ,D I  5. 5//// /  ) 
GC     TO     10CO 

400     STC° 
ESC 

SL3R0UTINE     E V AL ( A . F  .  AL  •  B  .T  ) 
IMPLICIT     RE*L»8     (A-h.O-Z) 
CI  OENSICN     i  (7  )  .»  (&  ) 
3  AT  A     Z/0.2C1194CS3  55  74  3  5D0  .  C. 39  4  IE  1 347  077563DO . 

1  0.57  09  7217  260e539DO,0 . 724 4 17 72 136 0  1  70D 0,  0. 8482065 834  I  04  2 7D C  . 

2  0.9  i 727 33 =24007  06  00.0 .987  99251 30  20  435  00/ 
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CAT*     V»/0  . 1534214 253271  HD0.C.l£61€lC0O1155e2D0« 

1  C.16£26;20;ei6S;4D0t0.1  39 57 C67 79261  j4DO,0.  l07  1592  20  467172COt 

2  0.07G366O47488108OO,0.O3O753241S9£U7D  Ot  C.  20257824  I  92556  1  C  0/ 
SUM=O.ODO 

C=(T-0  )/2. COC 
C=(T*E  )/2.0C0 
SCJ»=SCM+  »(?)«FUNCTN(D» 
DO  I      1=1.7 
Zl  =C*2( I  )*C 
ZZ=J-C*Zll ) 
.     SJM=SUM+W (  I  )* (FUNCTN(Z 1 )*FUNCTN(Z2) ) 
F=C«SUV 
HE  TURN 
ENC 

FUNCTION    FLNCTN(Y) 
IMPLICIT    REAL*6     (A-H.Q-Z) 
CCM^ON    X . AL 
V1  =  V(  I.  000  ,AL.Y) 
V2=V<  X.1.OC0.  Y) 
FUNCTN=(  1.0C0-V2/V  1  )**2 
RETURN 
ENC 


END  3F  FICE  CN  SVSIN. 
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THIS    PROGRAM     EVALUATES     THE     INTEGRALS     OF    SITE-SIYE     CIRECT 
CCRHELATCON     INTEGRALS     OF    THE    FORM     PROPOSED     BY    LOUDEN     AND     CHANDLER 
IMPLICIT     REAL*8     (A-H.O-Z) 

01  MENS  ION     ALOSOO)  .  RHOSOO)  .A  I  1  (  30)  .A21  (30)  ,A3l(30),A4l(30), 
*        A2  2 (30  )  ,  A32(30  ). A42(  30  ). A23( 30 J. A331 30)  tA43(  30 ) t A( 30) »T(  15) 
N        IS     THE    NLMBER    OF     STATE    CCNOITICNS    EVALUATED 
READ( 5.  1  )     N 

1  FORMAT (12) 

G*G      10     I     =      l.N 
10     R£AD(5. 2  )ALOS(  I >»RH03(  I  ) 

2  FCfiMAT(2Fl C.4) 
I     =      1 

*7     fiE^0{5,3)     A 

3  FORMA T(8F1 0.4 ) 

GO     TO     (  200.  20  1.  202,  2C3.2C4. 205. 20t .207,203.209)  .  I 
200     DC     100     J     =     l.N 
luC     A1KJ)     =    A(J) 

GC    T3     19 
20  1     3D     101      J    =     l.N 
10  1     \2  1(  J  )     =     A  (  J) 

GC     TO     19 
2  02    DC     1 C2     J     =     l.N 
10  2    A-3K  J  )    =    A(  J) 

GC     TO     19 
20  3    DO     I  03     J    =     l.N 
103     A4  I  (  J  )     =     A ( J  ) 

GC     TO     19 
204     OC      104     J     =      l.N 
10  4      A2  2  (  J  )      =     A  (  J  ) 

GO     TO     19 
20  5     DO     105     J    =      l.N 
1  J5     A32 (J )     =     A( J ) 

GC     TJ      15 
2J6     CO     106      J     =      l.N 

106  A42( J)      =     A ( J) 
GO     TO      19 

207     OC     107     J     =      l.N 

107  A23( J)      =     A ( J) 
GO    T]      IS 

203     CC     108     J     =      l.N 
106     A33( J )     =     A ( J) 

CO     TO     19 
209     OC     105      J     =     l.N 
109     A43( J )     =    A( J) 
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19     i      =     I     +     1 

1F(  I  .ZQ  .11)     GC     TO     2  1 

GC     TO    47 
<dl     JO     3  0     1  =  1  .N 

AL     =     AL3S(  I  ) 

3£     =     1.000     -    AL 

1(1)     =     (  A  I  1(  I  )-A21  (  D  +  A31  (  I)-A4  1  (  I  )  )/3.0D0 

T(2)      =      (  A21  (  I)-2.0D0*A3l  {  I  )*3.000*A41(  I  ))/4.0D  0 

T(3)      =     ( A3  1  (I )-3.0D0*A41  (I )  )/5.0DO 

T(  4  )     =     A4KD/C.CD0 

T (5)= (AL** 3)* (0.000    -AL*A22(I)+(AL**2> * A3 2( I )- ( AL* *3 > *A42 ( I > ) /3 .D 0 

T(6)=(AL**4)* (A22( I )-2  .  D0*AL*  A3  2  ( I ) +3. 00* ( AL**2) * A 42 ( I) >/4.D0 

T(7 )=(  AL  **5  )+( A32(  I  )-3. DO*AL*A 42 ( I)  1/5.00 

T(3)      =     ( AL**6 )*A42 ( I )/6.D0 

T(9)  =  BE**3MC.0DO    -<3E*A2  3(  I  )+ ( BE**2)*A33( I ) -( BE**3 > *A43  (  I  )  )/3. DO 

T  (10  >=BE**4*( A23(  I  )-2.D0*BE*A33(  I  )+3.D0*(  BE**2)*A43(  I  )  I/4.CC 

T( 11 )=BE**5*(A33( I  )-3 . DO*B E*A4 3 ( 1  ) )/5.00 

T(12)=(BE»*6)*A43(I )/6.00 

sux    =   o.oco 

DO     4  0     J     =      1.12 
40     SUM     =     SUM     ♦    T(J) 

C     =     16. 0D0*3. 141592*BH0S(I )*SUM 
AtO     =      l.OOC    -     C 

*fiIT£(6, 1000)C.AKO. ALOS( I).RHOS( I ).SUM 
1000    FOR^AM  1  OX,  •     THE    OCF     INTEGRAL     IS     '.FIG. 4.// 

*  10X,'         THE    DIMENSICNLESS     COMPREJS  IB1L ITV     IS        «.F15.4,// 

*  10X,«      L/S     =     '.F10.4,'         ANC       RhO*SlG**i    =     'iFIO.4// 

*  20X,      •         THE     INTEGRAL     iAS        '•     F20.6) 
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C  EXCXSS    ENTHALPY    CALCULATICNS    FRCM    LMFAC 

C  LIMITED     TO     FIVE    MOLECULES     »ITH     TEN    GROUPS 

DIMENSION        XM(5).X(10).R(10).a(lO).J»(l0.10).AOT(  10.   10) 

DIMENSION     AEXP< 10. 10) .ANUC 10.5) 

CIMENSION    SUMX(IO) 

COMMON     SUM 

R  =  AJ< 5. 123 )     NG.NM 
12.5     FCPMAT (211 ) 

REAlH  5. 124)  (C(I )  .1  =1  .NG) 

QC     100     1= l.NG 

READ(5. 1 24  )<A{I. J). J  =  l.NG> 

REAQ<  S, 124) ( ANU( I . J) ,J=l iNN) 
124     FORMAT(SFIO.O) 
100     CGNTINUE 

RSAD<5.*)     T 

CALL     ATEfiMS (A.NG.AOT,  A  EXP,  T) 

DC     1     1  =  1  .2  I 

AIM= 1-1 

XM1  ) =0.05*A  I* 

XM(2>  =1.  0-XM<  1) 

CALL     GFRAC (XM.X, NG.NM, ANU) 

CALL     hEG(X,G.ACT,AEXF,T .NG.HEXG) 

CALL     HEGR(  >«  .a.AOTi  AEXPiANU.T  iNN^C  tHEXKJ 

h£XGC=SUV»hEXC 

HE  *=HEXGC-HEXP 

*rite(6.  10)(XM(  I  I),  11=  l.NM) 
U     FCSMAT (1  OX,"  MCLE    FRACTIONS • 20X , 4 ( SX, F  10.4 )  ) 

»RITE( 6,  11 )  (  X(  I  I ),  I  1=1  ,NG) 

11  FORM  AT (  10X,  ■     GROUP     FRACTIONS •  ,  2  OX.  4(  5X,F  1  0.  4)  ) 

•RITE(6,12)     HEXGC.HEXfi.hEX 

12  FORMAT<  10X,  •     GROUP     CONTRIBUTION    =     '  .FlO.*i// 

♦  ,10X,      •         FEFERENCE    CONTRIBUTION        =         '.F10.4,//, 

*  ldX.'        E>C£SS    ENTHALPY        =        ".F10.5) 
1     CONT  INUE 

STCP 
END 

SUE  ROUT  INE     ATERMS( A.NG, AOT. AEXP. T ) 

JIM5NSI0N        XM(5).X(10).R{10).C(10).A(10.10).AOT(10.10) 
DIMENSION     AEXP(  1 0,  1 0). ANUl  10.5) 
JC1      1=1. NG 
DO     1     w=l .NG 
ACT(  I.  J )  =  A(  1, J )/T 
I     AEXPt  I, J)  =  EXF(-AOT  (  I, J  )  ) 
RETURN 
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SL3R0LTINE  GFRAC ( XM • X , NG ,SH  .ANU ) 

DIMENSION   XM(5).X(10).R(10).Q{  10) .A(  10,1  0) ,AOT(  10.10) 

JIMENSICN     AEXPdO  .10  )  ,  ANU{  10. 5i 

3IM£NSI0N     SLMX( 10) 

CCKMCN     SUM 

SL*=0.  0 

DC     2     1-l.NC 

SU*X( I )=0.0 

bo     1     J=1,NX 

5UMX(I)     =     S'JMX(I)      ♦     XM  (  J  )*ANU<  I.  J  ) 

SLM     =     SUM     +     SUHX (I ) 

QC     J     1=1. NG 

X( I )     =    SUMX  {  I  )/SUM 

RE  T  UR  N 

END 

SUBROUTINE     h£G ( X ,Q . A.T . N, HEX ) 

JIM£NSION     X(10).a(10).A(10.10).  AOT  (10.1  0).  AEXPdO.  10  ) 

30      I      1=1. N 

CC1J=1 .N 

AuT<I   .J)      =     A(I  . J)/T 

«eX?(I«J)     =     EXP(-AOT( I . J) ) 

SLI»A     =0.0 

00     2     K=1,N 

3UM3=0  .0 

SLMC=C«0 

CO     3    L=1.N 

ACC     =    X  (L  )  *C(L )*AEXP(L.K  ) 
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S0M3=SbM3+ADD*AC"I(L  ••<) 
3     SCMC=SUMC  +  ACD 

2  SLMA=SUMA+X ( K ) * Q { K ) *SU M8/SUMC 
t-.EX  =  S  JMA  *8.2  14*T 

RETJfiN 

ENiJ 

SLUR OUT INE     t-EGR(  XM,  Q. A OT, AEXP. ANU . T ,NM , NO »HEXR ) 

JICENSIGN        XM  (5) .X  C  10). R(l 0  ), 0(  10  ). A(  10. 10) .AOT(    10.10) 

DIMENSION     AE  XF  <  1  0,  1  0)  .  ANUU  0  .5  > 

SUMM=0  .0 

DC     5      1=1  iNN 

£CMA  =  C.  0 

DC    Z     K=l. NG 

SLMB=0.0 

SOMC  =  0. 0 

CC     3     L=l.NG 

ACC=ANU< U.  I  )*C(L)*AEXP (L.K) 

SUM3=SUMB+ADD*A0T(  L.K) 

3  5UXC=SUMC+/CC 

2  SLMA  =  SL'MA+AM,(K  .  I  )  *  C  (K  )  *SUMe/SU  MC 

3  SUMM=SUM*  +  SV.MA*XM(   I  ) 
HEXR==LMV*8.314*T 
HET  JRN 

cNC 

ENO     CF    F1L5     CN     SYSIN. 


APPENDIX  7 
COMPRESSIBILITY  THEOREM  FROM  RISM  THEORY 


The  purpose  of  this  appendix  is  to  show  that  a  general 
compressibility  theorem  exists  for  systems  for  which  the 
RISM  theory  applies. 

The  starting  point  for  this  derivation  is  the  constant 
temperature  Gibbs-Duhem  equation 

dPB  =  I   p .  dBy  .  (A7-1) 

i   x 

We  wish  to  expand  the  differentials  on  the  right-hand  side 
of  this  equation  as 

8Bu. 
dPB  =  11    p.  [j- L)         dp  (A7-2) 

ij  T,pWj     ^ 

And  then  we  must  identify  the  partial  derivatives  of  interest 
here.   We  do  know  that 

&  =px   [6i:+pJHi3] 

And  forming  these  elements  into  a  matrix  representation 
gives 
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And  then  a  thermodynamic  identity  gives 


-    Ufi   +    fiSfi]"1}^  (A7-5) 


>Pj   T,p 


Note  here  that 


!A7-6a 


(A7-6b: 


0.    +  £2p_  -  [I  +  £H]g 


.".[&  +  a  gg]  x  =  g  1  El  +  ggl  X 


and  the  ij-th  element  is  —  {[J  +  gg]   } — 
Thus 


dPB  -  H  (tl  +  all  1lii  dPj  (A7_7) 

ij 


So  now  we  must  determine  the  term  in  brackets 
To  do  this,  first  note  that 


[I  +  p_H]  .  .  =  6.  .  +  p  .  H.  .  (A7-! 


These  terms  must  be  related  to  the  site-site  correla- 
tion functions  of  the  RISM  theory.  Here  the  nomenclature 
of  Lombardero  and  Encisco  (1981)  is  used.   The  multicomponent 
RISM  equation  is  (in  Fourier  space) 
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h  *  wc  [I  -  PXc]  X  w  (A7-9) 


All  of  the  terms  here  can  be  considered  as  either  standard 
square  matrices  of  dimensions 


I  number  of  sites  in  molecule  i 

all 

types 

of 

molecules 

i 


or  simply  as  fourth  order  tensors.   To  explain  the 
nomenclature 


ha   =  total  correlation  function  between  site 
1^         a  on  a  molecule  of  type  i  and  site  8  on 
a  molecule  of  type  j . 


There  is  of  course  an  analogous  nomenclature  for  the  site- 
site  direct  correlation  functions. 
The  oo  matrix  has  elements  of  the  form 


)aB    =    6    ^aB  (A7-10) 

ij     i]   D 


where  co   is  the  Fourier  transform  of  the  intramolecular 


j  .  ; 
1 

correlation  function  for  sites  a  and  B  on  a  molecule  of 
type  j.   The  X  matrix  has  elements 
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vaf?  =  X.  6-  ■  ^6  (A7-11) 

xid    :  °i]  d 


We  shall  also  require  one  other  auxiliary  matrix,  a,  defined 


oaS  =  6. .  6  .  p .  (A7-12) 


Now  consider 


w  +  <?h].  .  =  m.  .  +  p,  h.  .  (A7-13) 


This  matrix  is  of  interest  because  (Chandler  and  Richardson, 
1984) 


k  iH    U  +  ohlj*  =  6±j  +  p±  Hi:j  (A7-14 


And  these  are  the  terms  required  to  construct  the  compres- 
sibility relation.   We  now  wish  to  relate  these  to  the 
site-site  direct  correlation  functions.   This  is  shown 
below. 

Using  the  RISM  equation  we  have 


u  +  gh  =  {I  +  oojC  [l-p$c]   }y  (A7-15) 

And  this  leads  to 
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I  +  Ohio  1  =  [1-pxc+atoc]  [i-pxc]  (A7-16) 


However,  a  simple  computation  will  show  that 

ou  =  px  (A7-17) 

And  thus 


+  oh  =  [I-pxc]  1    to  (A7-18 


Let   M  =  [I-pxc] 


Then  we  can  show 


Urn  r."  ,  _Ci<xB  _  y  C,a6,„,  =  Ma   _  ,Ma. 


:A7-19 


k,^S  {w  +  ah}l°    =    I    M^(o)  e  Mjj  =  (*£)..     (A7-20 


And  thus  we  have 


dPS  =  H    [(mV1].,  dP,  (A7-21) 

ij         13    D 

We  then  wish  to  determine  M   and  its  inverse.   From  the 
definition  of  the  matrix  inverse  we  know  that 


_1  (A7-22a 


I  =  MM 


or  in  component  notation 
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6.   .6  =    TV    MUY    (M"VB  (A7-22b) 

*Y 

But       M         =    I    -    PXC  (A7-23) 


Thus     (M"Vj    -    ,       «  -    p     J    Xl    SJ*    £« 


Then  we   have 

V.B    ■  ^  Mu  1  Vyb    -  p  1  ^  »?  S°8  '        <A7"25a) 

-  ^l  -  p  m  h  sjBt  »l6  ^  'A7-25b» 

And  in  the  long  wavelength  limit  this  is 

6.  .6  .   =  m"6(o)  -  p  HI  X„  M?f(o)  S"(o)  (A7-26) 

If  this  expression  is  the  summed  over  B  we  find 


).  .  =  M.  .  -  P  I   X   M.  .  II    C.  . 


(A7-27 


Define  C„,  =  II  C   (o)  (A7-28 


Then  in  matrix  form  this  is 


I  =  M   -  M  PC  lA/  " 
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a  -1 

And  therefore  M   =  (I-oc)  •         (A7-30 


and  then  (Ma)  1    =    I  -  pc  (A7-31) 


This  then  leads  to 


dP8  =  H    (6    -  p.  c..)dp.  (A7-32) 

ij 


And  we  then  have  the  compressibility  theorem 

(4-^1     =  11    X.  (6.  .-p.C.  .)  (A7-33) 

This  is  exactly  the  same  as  the  molecular  form  if  one  makes 
the  identification 


C  .  =  yy  CaB  (A7-34) 

1D   ae  13 


which  has  been  used  throughout  this  work.   It  should  be 
noted  here  that  this  derivation  requires  only  the  evaluation 
of  the  sum  of  site-site  functions  in  the  k+o  limit  which 
is  believed  to  always  be  well  bounded. 
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